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1  SUMMARY 

The  goal  of  the  Quantum  Computing  Approach  to  Model  Checking  for  Advanced  Manufacturing 
Problems  (QCHECK)  research  project  was  to  determine  if  it  is  feasible  in  the  future  to  speed  up 
a  Model  Checking  (MC)  approach  based  on  Counter-example  Guided  Abstraction  Refinement 
(CEGAR)  by  using  a  D-Wave  open  system,  adiabatic  quantum  annealing  processor.  These  are 
specialized  computing  devices  that  solve  spin  Ising  models,  which  are  equivalent  to  Quadratic 
Unconstrained  Binary  Optimization  (QUBO).  We  focused  on  two  aspects  of  the  CEGAR 
approach  that  involved  solving  integer  linear  programs  (ILPs)  and  Boolean  satisfiability  (SAT) 
problems.  The  project  was  divided  in  five  tasks: 

Taskl:  Study  feasibility  of  using  a  D-Wave  to  solve  Bounded  Model  Checking  (BMC)  problems 
and  implementing  Binary  Decision  Diagrams  (BDD)  based  techniques. 

Task  2:  Compare  the  performance  of  a  second  generation  D-Wave  (DW2)  on  MAX-2-SAT 
problems  native  to  its  architecture,  versus  the  heuristic  solver  MaxWalkSat. 

Task  3:  Develop  a  heuristic  embedding  algorithm  for  the  DW2  to  get  around  the  limited 
connectivity  of  the  processor. 

Task  4:  Integrate  the  CEGAR  approach  with  the  DW2  processor. 

Task  5:  Implement  examples. 

Results: 

Task  1:  It  was  found  that  even  though  in  principle  the  required  BMC  problems  could  be  cast  as 
QUBO  problems,  the  probabilistic  nature  of  the  processor  (that  provides  no  guarantees  that  the 
best  possible  solution  has  been  found)  made  the  approach  susceptible  to  false  negatives:  a  SAT 
formula  could  be  proclaimed  “unsatisfiable”  because  the  best  solution  found  by  DW2  does  not 
satisfy  the  formula,  while  a  better  solution  might  exist  that  proves  the  formula  satisfiable.  With 
regards  to  implementing  BDD  based  approaches  using  DW2,  it  was  concluded  that  encoding 
such  a  problem  as  an  optimization  problem,  though  possible,  would  not  scale  well  with  system 
size. 

Task  2:  The  comparison  was  performed  on  a  set  of  random  instances  of  MAX-2-SAT  that  are 
native  to  the  DW2  processor’s  architecture,  for  different  numbers  of  variables.  The  performance 
of  DW2  was  shown  to  be  better  than  that  of  MaxWalkSat,  with  the  caveat  that  MaxWalkSat  was 
not  optimized  for  the  DW2  architecture.  The  issue  of  DW2’s  performance  vs.  that  of  classical 
solvers  remains  (as  this  report  is  being  written)  an  open  and  very  contested  research  topic. 

Task  3:  A  tool  to  perform  heuristic  embeddings  was  created.  It  allowed  us  to  implement  QUBO 
problems  that  had  a  different  connectivity  graph  than  the  native  DW2  architecture.  This  tool  is 
useful  for  optimization  problems,  but  suffers  the  same  limitations  found  in  Task  1  for  decision 
problems  (i.e.,  SAT)  due  to  the  lack  of  a  guarantee  that  the  solution  found  is  the  best  possible. 
The  tool  needs  further  optimization,  and  alternative  approaches  need  to  be  investigated. 

Task  4:  Several  issues  arose  that  made  this  task  more  challenging  than  was  anticipated: 
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1 .  Off-the-shelf  model  ehecking  packages  do  not  provide  access  to  internal  data  such  as  the 
ILP  needed  to  be  solved  during  the  CEGAR  implementation.  We  contacted  the  developers  but 
they  were  reluctant  to  give  us  access  to  the  source  code. 

2.  Many  simple  models  that  we  studied  generated  trivial  ILPs  during  the  CEGAR 
implementation,  i.e.,  the  solution  could  be  found  by  simple  inspection.  It  took  some  time  to  learn 
from  these  models  what  properties  of  the  system  will  lead  to  non  trivial  ILPs. 

3.  In  order  to  have  access  to  the  ILPs  we  ended  up  writing  a  small  model  checking  package 
using  publicly  available  libraries.  We  used  the  And-Inverter  Graph  (AIG)  format,  allowing  us  to 
generate  small,  non-trivial  ILPs  (on  the  order  of  tens  of  linear  constraints  and  binary  variables). 

Task  5:  We  found  an  example  from  the  literature  of  a  flight  control  system  and  checked  a  safety 
property  that  requires  two  exclusive  flight  modes  not  to  be  engaged  at  the  same  time.  We 
implemented  this  problem  integrating  the  AIG-based  model  checker  with  the  DW2,  ran  the 
CEGAR  approach  starting  with  an  abstraction  that  had  33  hidden  variables.  Our  integrated  code 
proved  the  system  to  be  safe  by  making  visible  only  14  of  the  33  hidden  variables. 

Main  lessons  learned  during  the  execution  of  this  project 

The  current  programing  paradigm  of  the  DW2  processor  requires  either  heuristic  embeddings  or 
approximate  embeddings  to  implement  QUBO  problems  that  do  not  have  connectivity  native  to 
the  processor.  This  step  leads  to  a  loss  of  certainty  about  whether  the  optimal  solution  to  the 
original  problem  is  the  same  as  the  optimal  solution  to  the  embedded  problem.  This  feature 
potentially  results  in  false  negatives  when  solving  decision  problems  with  DW2.  Lor  example,  if 
the  best  answer  provided  by  the  processor  corresponds  to  a  negative  result  for  the  decision 
problem  (e.g.,  a  SAT  formula  is  not  satisfiable),  yet  there  exists  a  better  solution  to  the  original 
decision  problem  that  gives  a  positive  answer  (i.e.,  the  SAT  formula  may  be  satisfiable).  Note 
that  there  are  no  false  positives,  since  a  positive  answer  provides  that  the  corresponding 
assignment  can  be  checked  efficiently. 

Optimization  problems  are  better  suited  for  the  current  programming  toolbox.  Even  though  we 
may  not  find  the  optimal  solution  to  a  problem,  both  the  heuristic  and  approximate  embedding 
approaches  provide  “good  solutions”,  which  can  still  be  very  valuable  if  they  can  be  found  faster 
than  with  other  methods.  The  heuristic  embedding  tool  developed  in  this  project  is  designed  to 
generate  a  sequence  of  improving  solutions,  although  there  is  no  guarantee  that  the  optimal 
solution  will  be  found  (although  additional  information  about  the  problem  may  help  identify 
when  optimal  solutions  are  found). 

The  question  of  speed  up  with  respect  to  classical  algorithms  is  very  hard  to  answer.  A 
benchmarking  against  a  particular  classical  algorithm  will  not  preclude  the  existence  of  another, 
more  efficient  algorithm.  Since  we  can  only  estimate  the  scaling  behavior  of  the  runtime 
performance  of  the  DW2  processor  by  benchmarking  it  on  a  given  set  of  instances,  the  problem 
is  translated  into  finding  a  particular  set  of  instances  that  show  speedup  over  some  set  of  classical 
algorithms.  Even  how  to  pose  the  question  of  speedup  has  been  the  subject  of  intense  research. 
At  this  point  in  time,  there  is  no  conclusive  evidence  that  the  DW2  provides  any  speedup,  but 
this  has  only  been  tested  up  to  500  variables.  New  processors  with  up  to  2000  qubits  are 
expected  to  be  available  in  the  next  two  years. 

Although  not  directly  related  to  this  project,  very  important  results  have  been  obtained  regarding 
the  quantum  nature  of  the  DW2  processor.  Even  though  it  is  designed  to  operate  in  a  quantum 
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mechanical  regime,  it  is  not  easy  to  experimentally  confirm  this  feature.  During  the  execution  of 
this  project  we  also  performed  research  aimed  at  resolving  this  issue.  Two  approaches  were 
devised:  one  provided  evidence  of  a  quantum  signature  by  analyzing  the  statistics  of  the  output 
of  the  DW2  processor  when  solving  a  carefully  designed  problem  involving  8  qubits.  The  second 
approach,  implemented  in  collaboration  with  the  company  D-Wave,  gave  a  definitive  answer 
regarding  the  quantum  nature  of  the  device  by  showing  that  entanglement  is  present  during  the 
quantum  annealing  evolution.  Whether  this  entanglement  can  provide  a  computational  speedup  is 
still  an  open  question. 

In  terms  of  the  integration  of  the  CEGAR  model  checking  approach  and  the  DW2  processor,  the 
proposed  approach  was  shown  to  be  very  straightforward.  The  obstacles  encountered  were  not 
related  to  the  fundamental  idea  of  the  approach,  but  rather  to  the  technical  limitations  of  the 
software  tools  required  (lack  of  access  to  the  inner  workings  of  the  CEGAR  implementation 
available  in  the  different  publicly  available  model  checking  packages).  Any  model  checking 
package  that  provides  the  required  information  (i.e.,  the  ILPs  to  be  solved  in  CEGAR)  could  be 
easily  integrated  to  interface  with  the  DW2  processor. 
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2  INTRODUCTION 

The  goal  of  the  QCHECK  projeet  was  to  analyze  the  feasibility  of  exploiting  the  eomputational 
eapabilities  of  the  DW2  adiabatie  quantum  proeessor  in  order  to  speedup  and  improve  the 
solution  of  model  checking  problems.  The  DW2  device  is  designed  to  solve  combinatorial 
optimization  problems  by  exploiting  quantum  mechanical  effects  of  an  array  of  Superconducting 
Quantum  Interference  Devices  (SQUIDs)  [1]. 

One  of  the  main  drivers  of  the  computational  hardness  of  model  checking  problems  is  the 
extremely  large  size  of  the  state  space  that  needs  to  be  considered  [2],  The  different  algorithms 
and  techniques  that  have  been  developed  to  solve  model  checking  problems  need  to  implement 
in  one  way  or  another,  a  mitigation  strategy  for  this  problem.  One  of  the  approaches  that  have 
been  proposed  and  developed  is  based  on  abstractions.  The  main  idea  is  to  replace  the  system 
that  needs  to  be  checked  by  an  abstraction  that  has  a  much  smaller  state  space,  with  the  feature 
that  if  a  property  is  found  to  be  true  in  the  abstraction  it  is  automatically  true  in  the  original 
system.  Since  the  size  of  the  abstracted  state  space  is  smaller,  the  algorithms  employed  to 
address  the  abstract  problem  require  much  less  computational  resources. 

The  abstraction  based  approach  however,  comes  with  a  price;  a  property  may  be  proven  wrong  in 
the  abstraction  when  it  is  actually  true  in  the  original  system  (false  negative).  To  avoid  this 
problem,  every  counterexample  to  a  property  found  in  the  abstraction  must  be  verified  as  valid, 
i.e.,  a  corresponding  counterexample  must  exist  in  the  original  system.  When  such  a 
counterexample  cannot  be  produced,  we  say  that  the  abstraction  generated  a  spurious 
counterexample,  and  the  truth  or  falsehood  of  the  property  remains  unknown. 

To  solve  this  issue  an  approach  known  as  Counterexample  Guided  Abstraction  Refinement  has 
been  developed  [3].  The  basic  idea  is  to  use  the  structure  of  the  spurious  counterexample  to 
generate  a  finer  abstraction  that  would  get  rid  of  it.  A  finer  abstraction  has  a  larger  state  space 
and  so  it  is  important  to  find  a  refinement  that  increases  the  size  of  the  state  space  the  least.  This 
process  continues  until  the  property  is  proven  to  hold,  or  a  valid  counterexample  in  the  original 
system  is  found.  We  have  identified  an  approach  to  CEGAR  in  which  combinatorial  optimization 
problems  of  the  form  that  can  be  solved  by  the  DW2  processor  are  a  central  part  of  the 
algorithm;  one  is  to  check  whether  an  abstract  counterexample  corresponds  to  an  actual 
counterexample  in  the  original  system  which  requires  solving  an  instance  of  a  Boolean 
Satisfiability  problem;  the  other  is  at  the  root  of  finding  the  smallest  abstraction  refinement  that 
can  get  rid  of  a  spurious  counterexample  and  requires  solving  and  Integer  Einear  Program. 

2.1  D-Wave  Two  (DW2)  adiabatic  quantum  optimization  processor  overview 

The  DW2  adiabatic  quantum  computer  solves  a  Quadratic  Unconstrained  Binary  Optimization. 
This  optimization  consists  in  finding  the  vector  of  binary  variables  that  minimizes  the  quadratic 
objective  function 

fiXj,  ....  Xn)  ~  Min[^]  {  j=l,  }  (1) 

where  x  =  {xi,  X2,  ....  x„),  x,  e  {0,1},  and  Qijis  a  matrix  of  real  numbers  that  determines  the 
objective  function.  This  problem  is  equivalent  (through  a  simple  linear  transformation  of  the 
variables  Xi  Si  =  2xi  -  1)  to  the  Ising  model.  The  Ising  model  represents  a  set  of  interacting 
spin  magnets  with  an  energy  given  by; 
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E  {$1,  Sn)  i=l,  ...,n]hiSi  (2) 

where  the  spin  variables,  s'j’s,  now  take  the  values  {+1,-1},  the  parameters  Jy  represent  the 
interactions  between  two  spins,  and  the  parameters  hi  correspond  to  local  magnetic  fields. 
Solving  the  Ising  model  consists  in  finding  the  spin  configuration  that  minimizes  the  energy,  E. 
This  problem  is  known  to  be  NP-hard  [4],  and  many  important  combinatorial  problems  can  be 
reduced  to  it  [5]. 

DW2  implements  a  quantum  version  of  the  Ising  model,  where  each  spin  variable  is  replaced  by 
a  Pauli  operator  Cz,  representing  the  state  of  a  qubit  (quantum  bit)  that  is  associated  with  the 
magnetic  flux  of  a  superconducting  quantum  interference  device  (flux-SQUID).  The  Ising 
Hamiltonian,  given  by 

ffising  ~  CJj  ]}  +  (3) 

characterizes  the  quantum  mechanical  system  of  spins.  The  device  allows  for  tunable  interactions 
between  the  different  qubits  (i.e.,  tunable  parameters  Jy),  as  well  as  tunable  local  biases 
(parameters  hi). 

Quantum  annealing  in  the  D-Wave  processor  proceeds  as  follows:  initially  a  transverse  field  is 
applied  such  that  the  lowest  energy  state  has  all  the  spins  pointing  in  the  same  transverse 
direction,  a  quantum  superposition  of +1  and  -1.  The  parameters  are  then  slowly  varied  in  order 
to  transform  the  Hamiltonian  into  //ising,  whose  ground  state  encodes  the  solution  to  the 
optimization  problem.  The  adiabatic  theorem  of  quantum  mechanics  assures  us  that,  provided 
this  parameter  change  is  slow  enough,  the  final  state  of  the  system  corresponds  to  the  ground 
state  of  the  final  Hamiltonian  [6],  i.e.,  the  spin  configuration  that  minimizes  the  energy  function. 
The  values  of  the  spins  are  obtained  by  measuring  the  flux  of  each  qubit  at  the  end  of  the 
annealing.  In  reality,  due  to  the  probabilistic  nature  of  this  quantum  mechanical  system,  this 
process  must  be  repeated  several  times  in  order  to  expect  to  identify  the  lowest  energy 
configuration. 

Our  initial  D-Wave  One  (DWl)  system  had  128  qubits,  depicted  as  green  and  gray  circles  in 
Figure  1.  They  are  arranged  in  a  4  x  4  array  of  8-qubit  tiles.  In  each  tile,  the  qubits  are  separated 
in  two  groups  of  4  and  connected  in  a  bipartite  fashion  [7]  (each  qubit  is  only  connected  to  all  the 
qubits  in  the  other  group).  Some  qubits  in  each  tile  have  extra  connections  to  qubits  in  other 
tiles,  such  that  the  graph  is  connected  (but  not  fully  connected).  The  connectivity  graph  is  called 
the  Chimera  graph  [8].  The  DW2  processor  used  in  the  latter  part  of  the  project  has  512  qubits. 
It  is  composed  of  a  16  x  16  array  of  8-qubit  tiles  connected  in  a  similar  way  as  in  Figure  1 . 
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Figure  1,  D-Wave  One  connectivity  graph 


This  connection  topology  is  dictated  by  eonstraints  imposed  by  the  underlying  teehnology,  but 
the  design  is  sealable  up  to  many  thousands  of  qubits.  The  laek  of  full  eonnectivity  between  all 
the  qubits  in  the  ehip  prevents  a  straightforward  mapping  of  an  arbitrary  Ising  Hamiltonian  (or, 
equivalently,  an  arbitrary  quadratic  function)  into  the  proeessor.  However,  although  eonstructing 
and  optimizing  this  embedding  is  not  a  trivial  issue,  several  heuristies  have  already  been 
developed. 

2.1.1  The  Physical  Principles  of  the  D-Wave  Quantum  Computer 

The  basie  building  block  of  the  DW2  quantum  annealing  ehip  is  a  supereondueting  flux  qubit, 
(rf-SQUID  flux  qubit)  as  depieted  in  Figure  2.  The  simplified  version  eonsists  of  two 
supereondueting  loops  having  two  Josephson  junctions  [9].  Each  loop  is  subject  to  externally- 
biased  magnetic  fields  (Oix  and  02x)  that  are  used  to  control  the  properties  of  the  qubit.  The 
quantum  states  are  assoeiated  with  the  quantized  magnetie  flux  Oi. 
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Figure  2,  Schematic  representation  of  the  compound  superconducting  loops  used  to  realize 

the  quhits  in  the  D-Wave  processor 


For  low  temperatures,  it  is  a  good  approximation  to  only  eonsider  the  two  lowest  states 
eorresponding  to  flux  pointing  up  and  flux  pointing  down.  The  energy  profile  that  deseribes  the 
system  is  a  double-well  potential,  represented  in  Figure  3.  The  bias  fluxes  Oix  and  O2X  are  used 
to  adjust  the  height  of  the  energy  barrier  dU,  and  the  energy  difference  between  the  two  states, 
2h.  The  actual  qubits  in  the  quantum  computing  element  chip  have  extra  loops  that  are  used  to 
compensate  for  undesirable  effects  due  to  fabrication  manufacturing  variations  and  provide  more 
uniformity  in  their  properties. 


2.1.2  Programming  and  Using  the  D-Wave  Quantum  Computer 

Programming  of  the  DW2  involves  setting  the  values  of  the  local  magnetic  fields,  and  the 
coupling  coefficients  for  each  super-conducting  qubit,  which  determine  the  desired  final  (Ising) 
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Hamiltonian.  The  DW2  is  housed  in  an  integrated  environment  able  to  be  programmed  by  either 
desktop  computers  connected  directly  to  the  DW2,  or  remotely  via  a  local  area  network,  or  other 
remotely-accessed  communication  networks,  to  which  the  desktop  computers  are  connected. 

To  program  the  DW2,  a  user  needs  to  provide  the  chip  with  the  values  of  J  and  h  that  satisfy  the 
constraints  discussed  above.  Casting  a  given  problem  into  the  Ising  form  and  respecting  the 
constraints  on  J  and  h  is  part  of  the  “art”  of  programming  DW2.  For  many  discrete  optimization 
problems  there  are  known  mappings  to  the  Ising  model,  but  those  often  result  in  matrices  J  with 
more  connectivity  than  what  is  currently  available  on  DW2.  One  way  is  to  map  such  J  into  the 
chip  by  using  certain  qubits  to  simulate  more  connections,  but  the  price  paid  is  that  these  qubits 
are  not  available  to  encode  the  solution.  There  are  also  heuristic  approaches  that  aim  at 
approximating  a  given  unconstrained  J  with  another  matrix  J’  that  satisfies  the  constraints  of  the 
chip  and  has  the  same  minimum  of  the  energy  function. 

Programming  of,  and  readout  from  the  DW2  occurs  through  an  application-programming 
interface  consisting  of  function  libraries  that  make  calls  to  the  optimization  capability  of  the 
DW2.  These  libraries  are  available  in  Matlab  and  Python,  and  can  be  used  to  access  the 
machine’s  functionalities  directly  from  the  programmer’s  code.  These  software  tools,  in 
conjunction  with  the  machine’s  circuitry,  translate  the  description  of  the  Ising  Hamiltonian  into 
the  time-dependent  classical  controlling  signals  that  make  the  qubits  evolve  following  the 
required  adiabatic  evolution.  After  the  qubits  are  measured,  the  results  are  also  available  through 
a  software  interface. 


2.1.3  Counterexample  guided  abstraction  refinement 

One  of  the  main  computational  bottlenecks  in  model  checking  is  related  to  what  is  known  as  the 
“state  space  explosion”:  even  for  moderately  sized  systems,  the  state  space  needed  to  describe 
them  is  intractably  large  (a  system  with  10^°°  states  is  not  uncommon).  Developing  techniques  to 
deal  with  this  issue  is  central  in  model  checking.  One  popular  technique  is  based  on  the  use  of 
Binary  Decision  Diagrams,  a  very  compact  data  structure  that  allows  for  a  succinct  description  of 
the  state  space  and  the  transition  system  [10].  Another  approach  is  based  on  abstractions:  a 
smaller  system  is  constructed  in  such  a  way  that  properties  proven  true  in  the  abstraction  are 
guaranteed  to  be  also  true  in  the  original  system  [11].  The  abstraction  can  then  be  checked  using 
regular  model  checking  tools  (like  BDDs  for  example),  which  are  computationally  more  efficient 
since  they  are  applied  to  a  much  smaller  system. 

If  an  abstraction  is  not  sufficient  to  prove  a  given  formula,  the  model  checking  tool  used  on  the 
abstraction  must  provide  a  counter-example  (CE),  a  path  in  state-space  that  violates  the  formula. 
This  CE  can  be  real  or  spurious:  a  real  CE  can  be  mapped  to  an  actual  CE  in  the  original 
(concrete)  model,  hence  disproving  the  formula;  a  spurious  CE  is  an  artifact  of  the  abstraction 
and  “disappears”  when  mapped  to  the  original  model.  To  determine  which  one  is  the  case,  we 
can  “simulate”  the  CE  in  the  original  system.  This  can  be  posed  as  finding  a  satisfying 
assignment  for  a  Boolean  formula.  DW2  implements  these  Boolean  satisfiability  problems  by 
fabricating  an  energy  function  that  achieves  a  minimum,  when  all  clauses  are  satisfied.  If  the 
formula  is  not  satisfiable,  the  lowest  energy  configuration  obtained  with  DW2  will  represent  an 
assignment  of  the  Boolean  variables  that  will  not  satisfy  the  formula.  This  can  be  efficiently 
checked  from  DW2’s  output. 
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In  order  to  determine  if  an  abstract  CE  corresponds  to  an  actual  CE  in  the  original  system  (and 
hence  a  proof  that  the  property  is  not  satisfied),  we  need  to  translate  the  sequence  of  transitions 
that  form  the  abstract  CE  into  a  sequence  of  transitions  in  the  original  system.  The  central 
question  that  we  need  to  ask  is:  given  a  transition  in  the  abstract  system  does  a  corresponding 
transition  exist  in  the  original  systeml  It  should  be  remembered  that  in  an  existential  abstraction, 
an  abstract  transition  between  two  abstract  states  exists,  even  if  only  one  pair  of  original  states 
(each  mapped  to  a  different  abstract  state)  has  a  transition  (see  Eigure  4,  where  the  transition 
between  states  S3  and  Se  in  the  original  system  induces  a  transition  between  the  second  and  third 
states  of  the  abstraction). 

Eor  example,  consider  an  abstract  CE,  f,  given  by  a  sequence  of  abstract  states  {si , 

Given  an  abstract  state  s,  the  abstraction  function  h  maps  states  in  the  original  system  into  the 
abstract  system.  The  set  of  states  that  are  mapped  into  s  are  the  ones  that  satisfy  h{s)  =  s.  If  we 
denote  by  i?(Sj,Sy)  the  characteristic  function  of  the  transitions  in  the  original  system  (i.e., 
Sy)  =  1  if  and  only  if  there  is  a  transition  between  states  Sjand  Sj  ,  and  0  otherwise),  then  a 
path  (si, ... ,  Sji )  is  a  concrete  representation  of  the  abstract  CE  f,  if  the  Boolean  formula 

Af=i(/i(sj)  =  sJ  A  (4) 

is  satisfied,  where  the  first  AND  operator  assures  that  the  states  Sj  are  mapped  into  the 
corresponding  abstract  states  while  the  second  AND  operator  assures  that  there  exists  a 
transition  in  the  original  system  between  the  states  Sj  and  Sj+i.  We  will  show  later  that  finding  a 
satisfying  assignment  of  a  Boolean  formula  can  be  cast  as  a  0-1  Integer  Einear  Program,  and  this 
lEP  can  be  mapped  into  QUBO  form  implemented  by  DW2. 

2.1.4  Refining  the  Abstraction 

If  the  CE  is  spurious,  the  satisfiability  problem  presented  above  identifies  an  abstract  state  that 
causes  the  violation  of  the  formula  being  checked.  This  is  due  to  having  clustered  together 
''dead  end’’’  states  (that  do  not  provide  a  path  to  an  error  state)  and  "bad'  states  (that  provide  a 
path  to  an  error  state).  We  can  illustrate  this  behavior  with  the  following  diagrams  that  represent 
a  system  and  its  abstraction  (see  Eigure  4). 
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ORIGINAL  MODEL 


ABSTRACTED  MODEL 


Figure  4,  Example  of  abstracted  model 


Consider  the  set  of  states  and  their  transitions  depieted  in  the  original  model  of  Figure  4.  Here, 
So  is  the  initial  state,  and  the  states  eolored  red  are  error  states.  Consider  a  possible  abstraetion  of 
this  system  in  whieh  all  the  states  inside  each  dotted-lined  rectangle  are  mapped  to  the  same 
abstract  state.  The  transition  diagram  for  such  an  abstraction  is  given  in  the  lower  part  of  Figure 
4.  From  these  diagrams,  we  can  see  why  a  spurious  CE  may  arise.  In  the  original  system,  it  is 
clear  that  if  we  start  in  the  initial  state  I,  we  will  never  reach  the  “error”  states  (red  states  in  the 
diagram).  However,  by  looking  at  the  transition  diagram  for  the  abstraction,  we  can  see  that 
starting  in  the  abstract  initial  state  we  may  eventually  reach  an  “error”  state.  This  is  easier  to  see 
in  the  diagram  of  Figure  5. 


Figure  5,  Spurious  counter  example 
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The  states  in  black  in  Figure  5  represent  a  possible  path  that  the  system  can  take.  Note  that  in  the 
original  system  that  path  cannot  reach  any  “error”  states.  The  farthest  the  system  can  go  is  the 
“dead  end”  state.  However,  in  the  abstraction,  the  “dead  end”  state  is  mapped  to  the  same 
abstract  state  as  a  “bad”  state,  i.e.,  a  state  that  can  eventually  transition  into  an  “error”  state,  so  a 
spurious  transition  is  created  that  causes  the  system  to  have  a  valid  path  from  the  initial  state  to 
an  “error”  state.  In  the  diagram,  we  note  as  a  “failure”  state  the  abstract  state  that  is  causing  a 
spurious  counter-example  to  appear,  because  it  maps  together  “dead  end”  and  “bad”  states.  To 
refine  the  abstraction,  these  states  need  to  be  separated  and  assigned  to  different  abstract  states. 
Therefore,  in  order  to  construct  viable  abstractions,  one  must  follow  the  steps  outlined  below: 

1 .  Minimal  Separating  Set  -  Generate  an  abstraction  where  “dead  end”  and  “bad” 
states  are  clearly  separated.  This  results  in  an  ILP  that  the  DW2  can  solve  as  an 
Ising  problem. 

2.  Satisfiability  -  Regardless  of  the  approach  used  to  generate  abstractions,  any  CE 
needs  to  be  verified,  if  it  is  real  or  spurious.  This  invokes  the  processing  of 
satisfiability,  which  is  executed  in  DW2  as  an  Ising  model. 

The  sections  below  discuss  how  these  steps  are  formulated  and  mapped  to  the  DW2  quantum 
annealer. 

2.1.5  Minimal  Separating  Set 

This  can  be  constructed  by  identifying  which  of  the  “invisible”  variables  should  be  made 
“visible”  to  distinguish  “dead  end”  and  “bad”  states.  The  goal  is  to  separate  these  two  sets 
exactly  (i.e.,  no  mistakes  allowed)  using  the  smallest  number  of  invisible  variables  possible,  in 
order  to  keep  the  size  of  the  refined  abstraction  from  growing  too  much.  For  this,  we  use  the 
following  definitions,  assuming  that  there  are  2  sets  of  states  S  =  {si,  S2,  ....  Sm)  and  T=  {ti,  t2,  .... 
tn)  that  need  to  be  separated  {S  can  represent  the  “dead-end”  states  and  T  can  represent  the  “bad” 
states).  Let  fVhe  the  set  of  variables  required  to  specify  all  states  in  the  original  system. 

Definition-1:  A  set  of  variables  U  =  (uj,  U2,  .....  separates  S  from  T  if  for  each 

pair  of  states  {si,  tj),  Si^S  and  tj^T,  there  exists  a  variable  Ur  such  that  Si(ur)  7^  tj(ur). 

Definition-2:  Given  2  sets  S  and  T  per  Definition- 1,  find  the  smallest  set  of  variables  U  = 
{ui,  U2,  .....  Uk)^W,  that  separates  S  from  T.  The  set  U  is  called  the  minimal  separating  set. 

We  can  assign  a  binary  variable  v,-  to  each  variable  Ui ,  that  will  represent  whether  that  variable  is 
included  in  the  separating  set  or  not:  if  Ui  is  in  the  separating  set,  then  V;=l;  otherwise  is  zero  (the 
corresponding  variable  is  excluded).  The  conditions  are  that  for  each  {si,  tj)  pair  at  least  one  of 
the  variables  that  distinguish  between  the  two  states  must  be  selected.  Thus,  there  is  a  total  of 
mxn  conditions.  Under  this  formulation,  the  minimal  separating  set  can  be  solved  using  integer 
linear  programming  in  a  conventional  computer  to  attain  an  exact  solution. 

Minimize  Il|j=i,...,k]  vj  with  each  v,-  =  0  or  1  (5) 

Subject  to:  {\/seS)  (VteT)  V/  >  1  for  ^(Mj)  ^  t{u^) 

Assuming  that  S  represents  the  “dead  end”  states  and  T  the  “bad”  states,  the  objective  function 
aims  to  minimize  the  count  on  the  number  of  new  variables  we  are  including,  while  the 
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constraints  express  the  fact  that  every  pair  of  “dead  end”-“bad”  states  has  to  be  distinguished  by 
at  least  one  of  the  new  variables  we  are  ineluding. 

We  refer  now  to  the  case  that  the  sets  S  and  T  represent  the  “bad”  and  “dead-end”  states  denoted 
hy  Sb  and  So-  The  following  lemma  from  formal  methods  apply: 

Lemma  -  Let  Uhe  the  set  of  variables  separating  the  “dead-end”,  So,  and  “bad”,  Sb,  states.  Let 
an  abstraetion  function  h  ’  correspond  to  the  visible  set  F’  of  variables  realizing  the  abstraction. 
Also,  let  V  represent  the  entire  set  of  the  original  visible  variables;  then  the  following  applies:  F’ 
=  FU  U.  The  abstraction  function  h  ’  maps  So  and  Sb  onto  different  states  in  the  abstract  model. 

This  lemma  implies  the  following:  1)  the  number  of  visible  variables  have  increased  from  |  F|  to 
|F’|=|F|+|F|;  and  2)  using  this  augmented  set  of  visible  variables  the  abstraction  function  will 
map  Sb  &  Sd  to  different  abstract  states,  which  do  need  to  be  considered.  As  said  above,  the 
model  cheeking  tool  will  check  if  the  property  of  interest  is  valid  in  the  abstraction.  If  it  is  valid, 
then  the  procedure  is  over;  if  it  is  not  and  it  generates  a  spurious  CE,  the  process  needs  to  be 
repeated  until  a  set  of  visible  variables  is  identified  on  which  that  property  holds  or  a  valid  CE  is 
found. 
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3  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 


3.1  Feasibility  study  of  using  DW2  for  Bounded  Model  Checking  and  Binary  Decision 

Diagrams 

The  analysis  of  the  suitability  of  DW2  to  address  BMC  problems  and  to  implement  BDD-based 
model  eheeking  approaehes  was  a  theoretieal  exereise,  and  no  partieular  assumptions  were 
required.  We  will  present  the  details  of  the  problem  setup  together  with  the  results  of  this 
investigation  in  Section  4.1. 

3.2  Benchmarking  of  DW2  performance  on  MAX-2-SAT  against  classical  solver 

MaxWalkSAT 

The  main  goal  of  this  task  was  to  benchmark  the  performance  of  the  DW2  processor  on  a  native 
problem,  against  a  classical  heuristic  solver.  The  rational  for  this  task  was  that  benchmarking  up 
to  that  point  was  performed  against  classical  exact  solvers,  in  particular  AK-MAXSAT  [12]. 
Since  this  solver  provides  a  guarantee  of  optimality,  it  requires  more  resources  (i.e.,  run  time). 
The  DW2  processor  is  a  probabilistic  solver  (solution  is  provided  with  a  finite  probability  and  no 
optimality  guarantee  is  given),  and  so  it  is  not  fair  to  compare  it  with  exact  solvers. 

3.2.1  MAX-2-SAT 

The  choice  of  MAX-2-S AT  is  based  on  the  fact  that  this  problem  can  be  trivially  written  as  an 
Ising  problem  that  is  native  to  the  DW2  processor. 

Definition  of  MAX-2-SAT:  given  a  Boolean  formula  in  conjunctive  normal  form  with  2  literals 
per  clause,  find  the  maximum  number  of  clauses  that  can  be  simultaneously  satisfied. 

The  key  point  in  implementing  this  problem  using  the  DW2  processor  is  to  notice  that  for  each 
2-literal  clause,  we  can  construct  a  Hamiltonian  of  Ising  form  whose  ground  state  is  composed 
by  the  satisfying  assignments  for  the  clause.  For  example,  consider  the  2-literal  clause  (x,  v  xj). 
This  clause  is  satisfied  if  any  of  the  two  variables  is  TRUE.  To  map  this  problem  into  an  Ising 
form  we  will  associate  to  each  variable  the  state  of  a  qubit,  with  x  =  TRUE  |+1)  and  x  = 
EALSE  |-1).  Consider  then  the  2-qubit  Hamiltonian 

H  =  Oi  -  aj+  Oi  ®  ffj  )  (6) 

where  the  Gi  is  the  Pauli  operator  associated  with  qubit  i,  and  {|+1)  ,  |-1)}  are  its  corresponding 
eigenvectors.  Table  1  shows  the  truth  table  of  the  2-literal  clause,  and  the  energies  of  the 
associated  states. 


Table  1.  Truth  table  for  Boolean  OR  and  energies  of  associated  2-qubit  Hamiltonian 


Xi 

Xj 

Xi  V  Xj 

^7 

H 

E 

E 

E 

-1 

-1 

1 

E 

T 

T 

-1 

+  1 

0 

T 

E 

T 

+1 

-1 

0 
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T 

T 

T 

+1 

+1 

0 

We  can  see  that  all  satisfying  assignments  are  associated  with  states  of  energy  0,  while  the  only 
non-satisfying  assignment  corresponds  to  a  state  with  energy  1.  Hence,  the  ground  state  of  the 
Hamiltonian  H  is  composed  exactly  by  all  the  satisfying  assignments  of  the  2-literal  clause.  If 
any  of  the  variables  appear  negated  in  the  formula,  we  just  need  to  flip  the  sign  of  the 
corresponding  Pauli  operator  on  the  Hamiltonian. 

If  we  have  a  conjunction  of  many  clauses  Cja...aCm,  we  just  need  to  add  the  corresponding 
Hamiltonians  H  =  Hj  +  For  any  possible  truth  assignment  to  the  Boolean  variables,  the 

energy  of  the  state  will  be  increased  by  1  for  every  unsatisfied  clause.  Then  the  energy  of  the 
ground  state  of  //will  be  the  minimal  number  of  unsatisfied  clauses,  from  which  we  can  trivially 
infer  the  maximum  number  of  satisfied  clauses,  i.e.,  the  objective  of  the  MAX-2-SAT  problem.  It 
is  then  clear  that  we  can  look  at  MAX-2-SAT  as  a  native  problem  to  the  DW2  processor.  Even 
though  the  decision  problem  2-SAT  is  known  to  have  a  polynomial-time  solution,  the 
optimization  problem  MAX-2-SAT  is  NP-hard  (i.e.,  a  polynomial-time  algorithm  for  it  would 
imply  the  existence  of  a  polynomial-time  algorithm  for  all  problems  in  NP). 

3.2.2  MaxWalkSAT 

We  chose  the  MaxWalkSAT  [13]  solver  as  the  classical  algorithm  for  the  benchmarking.  This 
solver  applies  heuristic  methods  to  provide  an  approximate  solution  to  a  MAX-2-SAT  instance. 
Since  it  is  not  required  to  provide  any  guarantees  of  optimality  (as  exact  solvers  do)  it  can  run 
much  faster  on  many  instances.  We  considered  that  this  provided  a  better  comparison  between 
classical  solvers  and  the  DW2  processor. 

MaxWalkSAT  is  a  variant  of  WalkSAT,  a  heuristic  SAT  solver.  In  its  more  general  form, 
MaxWalkSAT  solves  the  weighted  SAT  problem,  in  which  each  clause  is  given  a  weight  and  the 
goal  is  to  maximize  the  total  weight  of  all  simultaneously  satisfied  clauses.  In  our  case,  we  set 
the  weights  to  I.  The  algorithm  for  WalkSAT  starts  with  a  random  truth  assignment  for  all  the 
variables,  then  randomly  selects  an  unsatisfied  clause  and  a  variable  within  that  clause  is  flipped. 
This  variable  can  be  chosen  either  at  random,  or  as  the  variable  whose  flipping  minimizes  the 
number  of  already  satisfied  clauses  becoming  unsatisfied.  So  in  a  sense,  it  is  a  mixture  of 
deterministic  local  search  and  random  jumps. 

3.2.3  Instance  ensemble 

In  order  to  implement  this  comparison  we  generated  MAX-2-SAT  instances  that  were  native  to 
the  architecture  of  the  DW2  processor.  We  generated  problems  with  N  variables,  for  N  = 
20,40, ...,500.  The  number  of  clauses  was  chosen  to  be  2N,  since  it  is  known  that  this  ratio  of 
clauses  to  variables  generates  instances  that  are  typically  hard  to  solve. 

The  ensemble  was  composed  of  1000  instances  for  each  value  of  N.  The  instances  were 
constructed  in  the  following  way.  For  each  value  of  N,  we  chose  N  qubits  that  formed  a 
connected  subset  of  the  processor  (to  avoid  assigning  variables  to  qubits  that  were  not  connected 
to  other  qubits  of  the  set).  Then  we  randomly  picked  M=2N  of  the  available  couplers  associated 
with  the  set  of  qubits  to  represent  the  2-literal  clauses.  Finally,  for  each  clause  we  randomly 
(probability  Vi)  negated  the  literals.  This  construction  assured  us  that  all  clauses  were  distinct, 
and  hence  the  total  number  of  clauses  was  indeed  2N. 
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3.2.4  Benchmarking  strategy 

To  compare  the  performance  of  the  classical  and  quantum  solvers  we  implemented  the  following 
strategy.  First,  by  using  the  exact  solver  AK-MAXSAT,  we  found  the  optimal  value  of  the 
objective  function  for  every  instance  in  the  ensemble.  This  value  was  later  used  to  estimate  the 
probability  of  success  of  each  solver  on  each  given  instance. 

Quantum  solver  (DW2):  we  ran  each  instance  a  thousand  times,  using  an  annealing  cycle  of  20 
microseconds.  We  compared  the  value  of  the  objective  obtained  in  each  run  with  the  known 
optimal  value,  and  used  this  information  to  compute  the  probability  of  success  for  each  given 
instance.  Then  we  used  this  information  to  compute  the  expected  number  of  repetitions  (or  runs) 
needed  to  obtain  the  optimal  value  at  least  once  with  at  least  99%  probability.  This  number  of 
repetitions  times  the  annealing  time  used  (20  microseconds)  was  the  performance  figure  we  used 
for  each  instance.  We  then  averaged  this  value  over  all  instances  with  the  same  number  of 
variables  N,  and  used  it  to  compare  with  the  classical  solver. 

Classical  solver  (MaxWalkSAT):  the  classical  algorithm  MaxWalkSAT  requires  another  input 
parameter  called  the  “cutoff’,  that  gives  an  upper  limit  on  the  number  of  iterations  performed 
before  stopping.  Clearly,  if  the  cutoff  is  small,  the  algorithm  will  be  faster  but  we  may  not  find 
the  optimal  solution.  On  the  other  hand  if  the  cutoff  is  very  large,  the  algorithm  will  take  more 
time  but  will  have  a  better  chance  of  finding  the  optimal  solution.  There  is  then  a  tradeoff 
between  the  value  of  the  cutoff  and  the  time  it  would  take  the  algorithm  to  find  the  optimal 
solution  with  probability  at  least  99%.  We  ran  every  instance  with  different  values  of  the  cutoff 
in  order  to  find  a  value  that  will  reach  the  optimal  solution  with  99%  probability  in  the  fastest  run 
time.  We  then  averaged  these  values  over  all  of  the  instances  with  the  same  number  of  variables 
N.  All  the  instances  were  run  on  a  Mac  Pro  with  a  2.6  GHZ  processor  and  48Gb  of  RAM. 

3,3  Development  of  heuristic  embedding  algorithm 

The  main  goal  of  this  task  was  to  develop  a  tool  that  would  allow  us  to  embed  problems  that  do 
not  match  the  processor’s  connectivity.  As  discussed  in  Section  2,  the  DW2  processor  has  a  very 
particular  connectivity  graph  called  the  Chimera  graph  that  is  the  result  of  design  compromises 
between  scalability  and  algorithmic  power.  The  connectivity  graph  is  sparse,  and  each  is  qubit 
connected  to  at  most  6  other  qubits. 

This  design  feature  has  an  impact  on  the  type  of  problems  that  can  be  embedded  in  the  processor. 
A  general  Ising  model  will  have  an  underlying  graph  of  couplings,  and  if  this  graph  is  not  a 
subset  of  the  Chimera  graph  we  need  to  implement  alternative  ways  of  embedding  the  problem. 
Even  if  a  given  instance  was  a  subgraph  of  the  Chimera  graph,  finding  the  appropriate  mapping 
is  an  instance  of  Subgraph  Isomorphism,  another  combinatorial  optimization  problem  that  may 
be  as  hard  as  the  original  Ising  instance.  Hence,  with  the  current  design  of  the  processor  we  have 
no  choice  but  to  develop  alternative  methods  to  embed  problems.  It  is  important  to  point  out  that 
this  issue  is  not  particular  to  the  D-Wave  processors.  For  any  implementation  of  adiabatic 
quantum  optimization,  the  connectivity  of  the  processor  will  be  associated  with  some  physical 
interaction  between  qubits.  These  interactions  tend  to  be  local  and  thus  require  that  the  qubits  are 
close  to  each  other.  This  will  put  a  limit  to  the  number  of  interactions  a  given  qubit  can  represent, 
since  the  number  of  local  neighbors  in  any  reasonable  architecture  will  be  limited  and  much 
smaller  than  the  total  number  of  qubits.  Hence  the  problem  of  embedding  is  central  to  the 
adiabatic  quantum  optimization  approach  and  not  just  a  D-Wave  issue. 
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Since  most  Ising  problems  of  interest  will  not  fit  direetly  into  the  Chimera  graph,  it  is  neeessary 
to  develop  teehniques  to  go  around  this  issue.  We  have  already  mentioned  that  an  exaet 
embedding  requires  solving  a  hard  problem  and  does  not  seem  to  be  a  sealable  solution.  Henee, 
we  need  to  apply  some  kind  of  approximate  embedding.  This  requires  disearding  some  of  the 
information  (i.e.,  the  eouplings)  that  defines  the  problem  in  order  to  generate  a  related  Ising 
problem  that  can  be  fit  into  the  Chimera  graph.  This  approaeh  has  been  used  by  Google  and  D- 
Wave  in  an  image  reeognition  application  [14],  where  they  redueed  the  training  of  a  strong 
elassifier  to  a  QUBO  problem.  Their  approximate  embedding  seheme  followed  a  greedy 
algorithm,  that  aimed  at  keeping  the  largest  eouplings  (in  absolute  value)  with  the  rational  that 
these  eouplers  will  be  more  important  in  determining  the  strueture  of  the  best  solutions.  It  is 
important  to  note  that  this  step  requires  a  eertain  preprocessing  of  the  input  instanee  that 
inereases  the  eomputational  resourees  required.  Also,  there  are  no  theoretieal  results  that  would 
guide  this  proeess  or  give  any  guarantees  on  the  quality  of  the  solutions  obtained. 

3.3.1  Iterative  heuristie  embedding 

In  order  to  address  the  drawbaeks  of  the  approximate  embedding  method,  we  eonsidered  a 
different  approaeh  that  aims  at  taking  advantage  of  the  sampling  eapabilities  of  the  DW2 
proeessor.  The  idea  eomes  from  an  approaeh  to  optimization  problems  known  as  “Probability 
oolleetives”  [15].  The  main  idea  is  to  replaee  an  optimization  problem  with  a  sampling  problem. 
Given  an  objeetive  funetion  over  binary  strings  G(x),  one  approaeh  to  find  its  minimum  will  be 
to  sample  from  its  Gibbs  distribution,  which  is  given  by 

P(x)  =  exp(-p*G(x))/Z  (7) 

where  beta  is  the  inverse  temperature,  and  Z  is  the  partition  funetion,  whieh  is  defined  as  Z  =  S 
exp(-P*G(x)),  where  the  sum  runs  over  all  binary  strings.  It  is  elear  that  this  distribution  is  biased 
towards  eonfigurations  that  have  a  small  value  of  G(x)  due  to  the  exponential  faetor.  So  if  we 
had  aecess  to  a  machine  or  algorithm  that  generated  samples  following  this  distribution,  with 
high  probability  we  would  obtain  the  minimum  eonfiguration. 

The  key  point  is  to  eonsider  the  DW2  proeessor  as  a  parameterized  sampler,  where  the 
parameters  are  the  loeal  fields  and  the  eouplings  {hiJij),  and  the  output  is  a  distribution  over  the 
set  of  binary  strings.  The  goal  is  then  to  find  a  set  of  parameters  that  produee  an  output 
distribution  that  is  “elose”  to  the  Gibbs  distribution  assoeiated  with  the  objeetive  funetion  G(x). 
This  sets  up  an  iterative  proeedure; 


1 .  Initialize  the  set  of  parameters  {hiJij). 

2.  Sample  the  output  of  the  DW2  proeessor  using  these  parameters. 

3.  Compute  measure  of  “eloseness”  between  this  output  distribution  and  P(x). 

4.  Update  the  parameters  {hiJij)  in  order  to  deerease  the  measure  of  “eloseness”. 

5.  Go  baek  to  Step  2. 

This  proeess  eontinues  until  a  termination  eriterion  is  reaehed.  Every  time  we  sample  the 
proeessor,  we  ean  eompute  the  value  of  the  objeetive  funetion  G  on  all  samples  and  keep  traek  of 
the  one  that  gives  us  the  minimum. 
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As  a  measure  of  “eloseness”  we  ehose  the  relative  entropy  (or  Kullbaek-Leibler  divergenee) 
between  the  two  distributions  [16].  This  measure  has  the  property  of  being  non-negative,  and 
vanishing  if  and  only  if  the  two  distributions  are  identieal.  If  we  eall  Q(x;  hiJij)  the  output 
distribution  of  the  quantum  proeessor  (that  depends  on  the  parameters  ),  the  relative 

entropy  between  Q  and  P  is  defined  by 

D(Q||P)  =  S  Q(x;  hJij)  log  (Q(x;  hJtj)  /  P(x))  (8) 

where  the  sum  is  taken  over  all  binary  strings.  Our  goal  then  is  to  find  the  values  of  the 
parameters  {hiJij)  that  minimize  the  relative  entropy.  This  is  an  optimization  problem  of  a 
eontinuous  funetion  over  a  set  of  eontinuous  variables,  so  we  ehose  a  gradient  deseent  method. 

In  order  to  eompute  the  eomponents  of  the  gradient  of  the  relative  entropy,  we  would  need  to 
know  the  funetional  form  of  Q(x;  hiJij).  However,  this  funetional  form  is  not  available  to  us,  and 
we  ean  only  sample  from  the  distribution  Q(x;  hiJij).  To  move  forward,  we  made  the  assumption 
that  Q(x;  hiJij)  was  the  Gibbs  distribution  assoeiated  with  the  Ising  energy,  that  is 

Q(x;  hiJij)  =  exp(  -P  Eising  (x;  hJij)  )  /  Zq  (9) 

where  ZQ{hi,Jij)  =  E  exp(-P*Eising  (x;  hiJij))  is  a  normalization  eonstant  that  depends  on  the 
parameters  QiiJij),  and  the  sum  is  over  all  binary  strings.  By  making  this  assumption  we  ean 
explieitly  eompute  the  eomponents  of  the  gradient  and  obtain 


Vj=  -p  {<(2x,  -l){2xj-l)  log(Q(x;  hJy)  /  exp(-p*G(x))))  -  <(2x,  ) 

( log(Q(x;  hJij)  /  exp(-p*G(x)))) } 


V,  =  -p  {<(2x,  -1)  log(Q(x;  hiJij)  /  exp(-P*G(x))))  -  <(2x,  -1))  <  log(Q(x;  hJy)  / 
exp(-P*G(x)))) }  (10) 


The  expeetation  values  that  appear  in  the  gradient  are  taken  with  respeet  to  the  distribution  exp(  - 
P  Eising  (x;  hiJij) )  /  Zq,  i.e.,  the  Gibbs  distribution  assoeiated  with  the  Ising  model  implemented 
on  the  proeessor.  Even  though  we  know  it’s  funetional  form,  this  expression  is  hard  to  eompute 
beeause  it  requires  summing  over  all  binary  strings  to  obtain  the  normalization  eonstant  Zq,  and 
this  sum  has  exponentially  many  terms,  making  it  impraetieal  for  large  problems.  In  order  to  get 
around  this  obstaele,  we  will  make  another  approximation  and  use  the  sample  averages  to 
eompute  the  expeeted  values.  The  sample  averages  ean  be  obtained  by  evaluating  the  expressions 
on  the  samples  produeed  by  the  proeessor.  Sinee  we  will  only  generate  a  fixed  number  of 
samples,  this  eomputation  can  be  done  efficiently. 

Note  that  the  algorithm  makes  two  approximations:  first,  it  assumes  that  the  output  distribution 
from  the  DW2  processor  is  a  Gibbs  distribution  in  order  to  compute  the  gradient  of  the  relative 
entropy;  and  second,  it  replaces  the  expected  values  over  this  Gibbs  distribution  by  the  sample 
averages. 
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3,4  Integration  of  CEGAR  approach  with  DW2 

The  CEGAR  algorithm  is  a  means  of  tackling  the  state  space  explosions  that  often  arise  in 
model-checking.  In  CEGAR,  one  initially  computes  an  abstraction  of  the  original  model  that 
can  be  model-checked  more  easily  than  the  full  model.  This  must  be  an  abstraction  that  is 
conservative,  in  a  sense  we  describe  below.  One  then  checks  the  abstracted  model  to  see  if  the 
property  holds  in  the  abstracted  model.  If  it  holds,  we  are  done;  the  system  passes  the  test.  Here 
is  where  the  conservative  nature  of  the  abstraction  is  critical:  it  must  be  the  case  that  if  the 
system  passes  the  check  the  property  is,  in  fact,  safe  (the  check  must  be  sound)',  however, 
CEGAR  admits  false  positives  (where  the  check  fails,  although  the  system  is  safe  -  the  check  is 
not  complete).  Typically  the  CEGAR  algorithm  is  applied  to  reachability  problems,  where  the 
safety  property  states  that  the  system  must  not  reach  some  undesirable  state.  A  conservative 
abstraction  is  used  which  increases  the  set  of  reachable  states,  so  that  the  check  will  be  sound. 

If  we  find  an  abstract  counterexample,  we  commence  the  part  of  the  process  that  gives  the 
algorithm  its  name.  Eirst  we  must  check  to  see  if  the  counterexample  is  sound.  We  do  this  by 
“replaying”  the  counterexample  in  the  full  model,  instead  of  the  abstraction.  If  the 
counterexample  is  found  to  be  sound,  we  are  done:  the  system  is  unsafe,  and  must  be  corrected. 
On  the  other  hand,  if  the  counterexample  is  unsound,  we  must  refine  the  abstraction  and  repeat 
the  process.  A  simple  diagram  of  the  CEGAR  procedure  is  presented  in  Eigure  6.  The 
abstraction  refinement  is  counter-example  guided  in  the  sense  that  we  find  a  place  in  the  counter¬ 
example  trace  where  the  abstract  counterexample  cannot  be  followed.  In  this  case,  what  must 
have  happened  is  that  the  abstract  counter-example  progresses  from  abstract  state  asi  to  asi+i  but 
there  is  no  way  to  progress  from  a  corresponding  concrete  state  h~^ (as i)  to  (as i+i)  (/z  '^is  the 
inverse  of  the  abstraction,  so  h~^(asi)  is  the  set  of  concrete  states  that  correspond  to  the  abstract 
state  asi.)  To  refine  the  abstract  state,  we  find  the  set  of  concrete  states  that  satisfy  the 
description  of  abstract  state  i,  h~^  (as  f  and  the  set  of  concrete  states  from  which  (asi+j)  is 
reachable,  and  refine  by  adding  state  features  that  separate  these  two  sets.  The  first  set  of  states 
is  called  the  “dead  end  states,”  and  the  second  set  is  called  the  “bad  states.” 


Figure  6,  CEGAR  loop 

In  the  QCHECK  project,  we  experimented  with  applying  the  D-wave  quantum  computer  to  two 
steps  of  the  CEGAR  process: 
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1.  Checking  the  soundness  of  the  abstract  counter-example  and 

2.  Finding  an  (approximately)  optimal  refinement  to  separate  the  dead-end  from 
the  bad  states. 


The  first  step  involved  solving  a  SAT  problem  that  encoded  the  existence  of  an  actual  counter 
example  associated  with  the  abstract  counterexample.  The  second  involved  solving  an  ILP  in 
order  to  find  the  smallest  set  of  hidden  variables  that  needed  to  be  made  visible  in  order  to  get  rid 
of  the  spurious  counterexample.  These  two  problems  (SAT  and  ILP)  are  special  cases  of 
combinatorial  optimization  problems,  and  can  both  be  cast  as  QUBO  problems  that  can  be 
solved  with  DW2.  Actually,  SAT  can  be  cast  as  a  particular  instance  of  ILP  as  we  show  below. 

3.4.1  Converting  SAT  to  ILP 

In  a  Boolean  satisfiability  problem,  we  have  a  set  of  Boolean  variables  and  a  set  of  clauses 
formed  by  combining  a  number  of  those  variables  and  their  negations  with  the  logic  operator 
OR.  For  example,  a  clause  can  take  the  form 

C  =  V  :t2  V  (11) 

where  V  stands  for  the  OR  operator  and  ^  stands  for  the  negation  of  :t:3 .  A  clause  is  satisfied  if 
at  least  one  of  its  Boolean  operands  is  TRUE.  The  Boolean  satisfiability  problem  consists  in 
finding  an  assignment  of  truth-values  for  all  the  variables  such  that  all  clauses  are  satisfied.  As 
we  have  shown  above,  a  step  in  the  abstraction  refinement  process  requires  us  to  solve  a  Boolean 
satisfiability  problem  in  order  to  determine  if  an  abstract  CE  corresponds  to  a  CE  in  the  concrete 
system.  We  will  now  show  how  to  map  this  problem  into  an  lEP.  Eor  each  clause,  we  will 
replace  the  OR  operator  by  the  “+”  (sum)  operator,  every  Boolean  variable  Xi  by  a  binary 
variable  x,  if  it  is  not  negated,  and  by  (1-  x,)  if  it  is  negated,  and  we  will  impose  the  condition  that 
this  sum  must  be  greater  or  equal  to  1.  We  will  use  the  numerical  value  1  to  represent  the 
Boolean  value  TRUE,  and  0  to  represent  EALSE.  Hence,  the  clause  C  defined  above  will  be 
transformed  into 

^1 +^2  +  (1 -^3)  ^  1  (12) 

It  is  not  difficult  to  see  that  the  logical  clause  C  is  TRUE,  if  and  only  if  the  above  linear 
inequality  is  satisfied.  Then,  satisfying  a  set  of  clauses  is  equivalent  to  checking  the  feasibility  of 
satisfying  a  set  of  linear  inequalities  like  the  one  above  (one  for  each  clause).  Since  the  variables 
are  binary,  this  is  nothing  but  a  particular  case  of  an  Integer  Einear  Program,  one  in  which  there 
is  no  linear  function  to  optimize  and  the  variables  are  restricted  to  be  0  or  1 .  This  type  of  problem 
can  be  implemented  on  DW2  as  explained  in  the  next  section. 

3.4.2  Converting  lEP  to  QUBO 

lEPs  are  known  to  be  computationally  expensive  for  conventional  computers  (since  they  are  NP- 
hard).  Nevertheless,  we  can  map  this  problem  to  an  Ising  model  that  can  be  implemented 
natively  on  the  DW2  quantum  processor.  Eet  us  show  in  general  how  this  mapping  can  be  done. 
A  typical  binary  integer  linear  programming  problem  can  be  written  as: 

minimize[x]  L|j=i,...,n]CjXj 

subject  to:  E|j=i^,,^n]5kAj  >  ;  hk>0  ;1<A:<^  (13) 
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where  eaeh  xj  is  a  Boolean  number  taking  the  value  of  0  or  1,  xje  {0,1},  all  the  eoeffieients,  ej, 
Bkj,  are  integer  numbers  (positive  or  negative),  and  hk  is  non-negative  integer  number.  The  way 
to  deal  with  these  constraints  is  to  transform  them  into  a  quadratic  penalty  term  that  will  increase 
the  energy  of  the  Ising  model  when  they  are  not  satisfied.  This  is  easily  done  with  equality 
constraints.  For  the  inequality  constraints,  we  first  transform  them  into  equality  constraints  by 
adding  extra  '"slack"  variables  (this  is  a  canonical  way  of  dealing  with  inequalities  in 
optimization  problems).  The  detail  we  have  to  take  into  account  is  that  the  variables  that  we  add 
must  be  binary  too  (as  the  xj’s).  For  the  k‘^  inequality  constraint,  we  define  mk  =  2[i=i,  ...,n]5ki 
including  only  the  5ki’s  that  are  non-negative  (5ki  >0).  In  essence,  mk  represents  the  maximum 
numerical  value  of  the  sum,  since  each  xj  is  a  Boolean  number.  The  role  of  the  slack  variables  is 
to  add  whatever  value  is  needed  in  order  to  transform  the  inequality  constraint  into  an  equality 
constraint.  The  slack  variables  must  construct  a  non-negative  integer  number  such  that 

2^D'=i,--,nAjXj-;?k  =  hk  (14) 

In  order  to  satisfy  this  equality,  the  integer  number  /»k  should  be  allowed  to  be  as  large  as  the 
maximum  difference  between  L|j=k...,  nj^kjxj  and  by:,  this  constitutes  its  numerical  range.  This 
maximum  difference  is  just  \my  -  bk\.  The  most  efficient  way  (in  terms  of  using  the  smallest 
number  of  extra  slack  variables)  is  to  express  the  non-negative  integer  py  by  an  expansion  in 
powers  of  2  (binary  expansion).  Therefore, 

=  (15) 

where  each  Pj  is  a  binary  variable,  Pj  e{0.1},  and  Dk  =  [|log2(mk-hk)|]+l,  since  py  <  \{my  -by)\. 
This  is  a  consequence  of  the  binary  arithmetic  stating  that  the  number  of  bits  required  to  express 
a  positive  number  of  magnitude  <  N  is  equal  to  log2N  +  1 .  Here,  the  number  of  bits  required 
(and  hence,  the  number  of  binary  slack  variables)  will  be  Dk,  and  the  square  brackets,  [  ...  ], 
represent  the  largest  integer  smaller  than  the  argument.  So  we  define  new  slack  binary  variables 
Pki  that  transform  the  inequality  constraints  in  Equations  (13)  to  equality,  namely, 

^|j=i,...,n]-SkjXj  -  5^[j=i,  ...,Dk]2-’  ^  /3j  =  by;  \  <  k<  q  (16) 

It  should  be  clear  that  if  this  equality  is  satisfied.  Equation  (16),  then,  the  inequality  in  Equation 
(13)  is  also  satisfied.  To  generate  an  Ising  problem  we  add,  to  the  linear  objective  function  of  the 
binary  ILP,  Equation  (13),  the  square  of  the  difference  of  the  left-hand-side  and  the  right-hand- 
side  of  each  equality.  Equation  (16),  times  a  penalty  constant  K>0,  and  we  now  take  the 
minimum  over  the  original  variables  Xi  and  the  new  slack  variables  pkj-  We  finally  get: 

njCjXj +K2j[k=l,...,  q]((2[j=l,...,  n].®kjyj  ■  ^[j=l,  Dk]2'^  flkj)  ~  b^  }  (17) 

where  x=  {xi;  i=,..,n},^=  {y^kj;  j=l,..,  Dk;  k=l,..,  q},  and  Dk  =  [|log2(mk-hk)|]+l  fork=l,..,q. 

Equation  (17)  has  the  form  of  an  Ising  model  (a  linear  term  and  a  quadratic  term).  We  can  make 
this  more  explicit  by  expanding  this  expression  and  grouping  the  corresponding  terms.  In  order 
to  make  things  more  compact  we  define  the  following  new  integer  matrices: 


-  ^ik 

if  1  <  i  <  n 

and  1  <  k  < 

q- 

d^ik 

_  2i-(n  +  X[l=  1,  ...,k-l]D/) 

if(n+  l<i< 

(n  +  E[i=i, 

k-i]D/)  and  1<  k<q. 

^=0 

Otherwise 

(18) 
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and 


=  Ci  if  1  <  i  <  n 

=  0  if(n+l<i<  (n  +  E[k=i, ...,q]Dyt)  (19) 

where  O  =  {0^  k=l,  q;  i=l,  n  +  i:[k=i,  is  a  qx  (n  +  i:[k=i,  ...,q]D;t)  matrix,  andg  is 

a  lx(n  +  2[k=i,  ...,q]D,t)  row  vector.  Also,  we  augment  the  original  vector  x  to  a  new  vector  z,  by 
appending  the  slack  variables  namely  creating  a  new  (n  +  2[k=i,  ...,q]D,t)xl  column  vector  ,  z,  z 
=  {xi,...,Xn,  pkj  ;k  =  l,...,q;  j  =  1,  ...,  Dk}. 

Equation  (17)  can  be  transformed  into  an  Ising  model  solvable  by  the  DW2  quantum  computer. 
The  minimization  is  over  the  augmented  vector,  z: 

min[,j  {  z^(0^0)z  +  (g  -  2b^O)z  }  (20) 

where  b  is  the  qxl  column  vector  {b\,  ...,  hq}. 

In  summary,  this  ILP -based  approach  to  abstraction  refinement  can  be  cast  as  an  Ising  problem 
that  can  be  natively  implemented  on  the  D-Wave  quantum  computing  processor. 

3.5  Implementation  of  Model  Checking  example 

The  final  task  of  this  project  consisted  of  applying  the  tools  and  techniques  developed  in  the 
previous  tasks  to  a  simple  model  checking  example.  We  considered  different  systems  and  finally 
converged  on  a  model  that  was  inspired  by  an  avionics  example  problem,  but  did  not  correspond 
to  any  real  hardware  or  software.  It  was  a  toy  model  designed  to  show  a  proof  of  concept  for  the 
integration  of  DW2  and  the  CEGAR  approach  to  model  checking,  and  was  tweaked  to  have 
certain  features  that  would  result  in  non  trivial  problems  for  the  DW2  processor  to  solve.  A  more 
detailed  description  of  the  example  will  be  provided  when  discussing  the  result  of  the  actual 
implementation  in  the  following  section. 
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4  RESULTS  AND  DISCUSSION 


4.1  Feasibility  study  of  using  DW2  for  Bounded  Model  Checking  and  Binary  Decision 

Diagrams 

Bounded  Model  Cheeking  is  an  alternative  approaeh  to  model  eheeking.  Its  main  feature  is  that 
properties  are  eheeked  to  hold  for  a  finite  number  of  time  steps.  In  this  way,  the  properties  ean  be 
expressed  as  Boolean  formulas  and  the  algorithm  eonsists  in  applying  a  satisfiability  solver  to 
determine  if  the  formula  has  a  satisfying  assignment:  if  it  does,  the  property  is  proven  true,  and  if 
it  does  not,  it  is  proven  false. 

This  approaeh  is  different  from  the  more  general  model  eheeking  approaeh,  whose  algorithms 
are  based  on  eonstrueting  a  sequenee  of  sets  that  tend  to  the  set  of  all  states  that  satisfy  a 
partieular  temporal  logie  formula.  The  advantage  of  the  general  approaeh  is  that  properties  ean 
be  proven  to  hold  for  all  possible  exeeutions  of  the  system,  while  BMC  are  only  restrieted  to 
exeeutions  within  a  finite  temporal  horizon.  The  priee  paid  in  BMC  is  a  laek  of  eompleteness, 

i.e.,  it  is  not  possible  to  prove  or  disprove  every  formula  in  a  given  temporal  logie.  However, 
there  are  eertain  properties  that  ean  be  proved  and  others  that  ean  be  disproved.  In  partieular, 
BMC  is  well  suited  for  finding  short  eounterexamples,  so  its  goal  leans  more  towards  finding 
bugs  than  proving  eorreetness.  In  this  area,  BMC  ean  be  more  effieient  than  general  model 
eheeking  teehniques  based  on  BDDs. 

BMC  proeeeds  in  two  steps:  first,  a  finite  length  execution  path  satisfying  a  certain  property  on 
the  space  state  is  encoded  as  a  propositional  formula;  then,  a  satisfiability  solver  is  applied  to 
find  a  satisfying  assignment  or  prove  none  exists.  If  a  satisfying  assignment  is  found,  it  can  be 
decoded  to  represent  a  particular  path  on  the  state  space  that  satisfies  the  property.  Depending  on 
the  property  being  considered,  this  could  be  a  proof  of  a  liveness  property  (i.e.,  a  state  with  a 
certain  property  can  actually  be  reached),  or  a  counterexample  that  disproves  a  property  (by 
showing  a  specific  path  that  violates  it).  It  is  in  this  second  step  of  BMC  that  we  believe  DW2 
can  provide  an  advantage,  since  Boolean  satisfiability  is  a  decision  problem  that  can  be  cast  as 
the  type  of  combinatorial  optimization  problem  that  DW2  is  designed  to  solve. 

4.1.1  Creating  propositional  formulas  in  BMC 
In  BMC  we  consider  three  elements: 


1 .  A  transition  system  M. 

2.  A  temporal  logic  formula  O. 

3.  A  time  bound  k. 


From  these  three  elements  we  construct  a  propositional  formula  that  checks  the  satisfiability  of 
the  property  represented  by  O  in  the  transition  system  M  for  paths  of  length  at  most  k.  We  can 
construct  the  unrolled  transition  relation  defined  by 

/(So)A  (21) 

where  /(sq)  is  the  characteristic  function  of  the  set  of  initial  states,  and  r(Sj,  Sj+i)  is  the 
characteristic  function  of  the  transition  relation.  Basically,  [M]^.  is  the  set  of  all  allowed  paths  of 
APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 

22 


length  k  that  start  in  the  set  of  initial  states.  We  can  also  construct  [0]/^,  which  is  a  formula  that 
will  be  true  if  and  only  if  O  is  valid  along  a  path  of  length  k.  In  BMC,  we  want  to  find  whether 
the  conjunction  formula  A  [0]/^)  has  a  satisfying  assignment. 

As  an  example,  consider  the  Computation  Tree  Logic  (CTL)  formula  EF  p,  that  means  that  a 
state  satisfying  the  propositional  formula  p  is  reachable  from  the  initial  state.  Applying  the  BMC 
approach  to  this  formula  for  k=2,  results  in  the  following  formula: 

/(so)  A  rCso.Si)  A  7’(Si,S2)  a  (p(So)  V  p(si)  V  p(s2)  )  (22) 

where  [EF  p]2  =  (p(So)  V  p(si)  V  p(s2)  ).  Solving  this  formula  consists  in  finding  a  sequence 
of  three  states  (Sq,  S2)  that  satisfies  it,  or  proving  that  no  such  satisfying  assignment  exists. 


4.1.2  Mapping  of  propositional  formulas  for  BMC  into  DW2 

Our  approach  to  implementing  a  satisfiability  solver  using  DW2  is  based  on  a  transformation  of 
SAT  to  a  0-1  ILP,  which  in  turn  can  be  transformed  into  a  QUBO  problem.  The  transformation 
from  SAT  to  0-1  ILP  assumes  that  the  SAT  formula  is  given  in  conjunctive  normal  form  (CNF), 
that  is 


F=CiA-AC^  (23) 

where  each  clause  Cj  is  a  disjunction  of  literals.  The  formulas  generated  in  the  BMC  approach 
need  not  be  of  this  particular  form.  Even  though  it  is  a  conjunction  of  three  elements  (the 
characteristic  functions  of  the  set  of  initial  states  and  the  unrolled  transition,  together  with  the 
formula  specifying  the  particular  property  being  considered)  each  of  these  terms  will  not 
necessarily  be  just  a  disjunction  of  literals. 

It  will  then  be  necessary  to  transform  them  into  CNF  before  we  can  map  the  problem  into  the 
chip.  This  type  of  preprocessing  is  also  present  in  classical  implementations  of  BMC,  and 
researchers  have  developed  subroutines  that  perform  this  translation.  The  drawback  is  that  they 
usually  require  the  addition  of  extra  variables,  which  will  reduce  the  number  of  qubits  available 
to  represent  state  variables,  resulting  in  a  reduction  of  the  size  of  the  systems  we  will  be  able  to 
analyze.  Given  that  the  number  of  qubits  available  is  fixed,  and  even  though  it  is  expected  that 
this  technology  will  scale  with  time  it  will  only  do  it  moderately,  it  is  clear  that  preprocessing 
techniques  and  hybrid  approaches  that  allows  us  to  break  the  problem  into  smaller  ones  will  be 
an  extremely  important  component  of  any  application  of  quantum  computing  in  the  adiabatic 
model  implemented  by  D-Wave  devices. 

4. 1 .3  Results  of  the  suitability  study  for  Bounded  Model  Checking 

The  analysis  presented  above  shows  that  the  BMC  approach  is  indeed  suitable  for 
implementation  using  the  capabilities  of  DW2.  Since  the  problem  reduces  to  SAT  and  SAT  can 
be  cast  as  a  combinatorial  optimization  problem  of  the  form  solved  by  DW2,  the  main  issues  to 
address  are  how  to  efficiently  use  the  resources  of  DW2,  namely  how  to  encode  the  SAT  formula 
in  a  way  that  it  makes  the  most  efficient  use  of  the  qubits  available.  These  are  not  different  from 
the  ones  that  we  had  already  anticipated  when  we  proposed  to  use  DW2  to  check  if  the  abstract 
counter  examples  obtained  within  the  CEGAR  framework  are  associated  with  real  counter 
examples  in  the  original  system  (this  is  again  an  instance  of  SAT). 
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4.1.4  Implementation  of  Binary  Decision  Diagrams  using  DW2 

The  purpose  of  this  analysis  was  to  study  whether  techniques  based  on  BDDs  used  in  model 
checking  of  finite  discrete  systems  can  be  cast  as  an  optimization  problem  that  would  be  suited 
for  implementation  with  DW2.  We  will  first  discuss  the  basic  features  of  BDDs,  why  they  are 
useful  in  MC,  and  the  issues  that  arise  when  we  try  to  implement  BDD  based  algorithms  as 
combinatorial  optimization  problems. 

4. 1 .5  Binary  Decision  Diagrams 

BDDs  are  a  type  of  data  structure  used  to  represent  Boolean  functions.  A  Boolean  function  is  just 
a  function  from  {0,1}^  to  {0,1},  that  is,  a  function  from  the  set  of  all  N-bit  strings  to  {0,1}.  We 
can  interpret  these  functions  as  the  characteristic  function  of  a  given  subset  to  N-bit  strings  (the 
subset  that  evaluates  to  1).  Since  the  model  checking  problem  is  essentially  a  set  problem  (i.e., 
whether  the  set  of  initial  states  I  is  included  in  the  set  of  states  that  satisfies  a  formula  O),  BDDs 
are  a  natural  tool  for  model  checking  algorithms. 

Why  are  BDDs  so  useful  in  MC?  The  main  problem  of  MC  is  the  state  space  explosion.  BDDs 
provide  a  structure  that  in  some  cases  can  represent  both  the  subsets  of  the  state  space  and  the 
state  transition  relation  in  a  compact  way.  Furthermore,  we  can  define  operations  on  BDDs  that 
implement  the  set  operations  required  by  the  model  checking  algorithms  (unions,  intersections, 
etc.)  This  allows  the  algorithms  to  operate  not  with  the  subsets,  but  with  their  characteristic 
functions.  The  point  is  that  a  subset  can  have  an  exponential  size  (in  the  number  of  variables 
used  to  describe  it)  while  its  characteristic  function  can  be  described  much  more  compactly.  In 
this  way,  the  algorithms  proceed  in  a  way  that  it  does  not  require  an  explicit  description  of  an 
exponentially  large  set  at  any  point. 

Formally  speaking,  a  BDD  is  an  acyclical  directed  graph,  with  one  root  node,  two  terminal 
nodes,  and  a  set  of  internal  nodes  that  have  one  predecessor  and  two  successors.  Each  node 
represents  a  variable,  and  the  value  of  the  Boolean  function  the  BDD  represents  is  obtained  by 
traversing  the  graph  starting  from  the  root  node,  and  moving  to  the  successor  node  associated 
with  the  value  of  the  variable  the  node  represents;  the  terminal  node  that  is  reached  is  the  value 
of  the  function  for  that  particular  bit  string. 

One  of  the  useful  features  of  a  BDD  is  that,  once  a  variable  ordering  is  fixed,  the  BDD  can  be 
written  in  a  unique  canonical  way,  which  makes  comparing  two  BDDs  (i.e.,  two  subsets  of  N-bit 
strings)  computationally  simple  for  polynomially  sized  BDDs.  The  size  of  this  canonical  BDD 
depends  crucially  on  the  variable  ordering,  and  two  different  variable  orderings  can  result  in  two 
BDDs  of  exponentially  different  size.  Even  though  finding  the  best  variable  ordering  for  a  BDD 
is  itself  a  hard  computational  problem,  heuristics  have  been  developed  that  result  in  reasonable 
sized  BDDs  for  many  problems  studied  in  practice. 

Another  important  reason  for  working  with  canonical  BDDs  is  that  the  set  operations  required  by 
model  checking  algorithms  can  be  implemented  directly  on  the  BDD  by  applying  a  small  set  of 
graph  operations:  given  two  BDDs  in  canonical  form,  we  can  compute  the  canonical  BDD 
corresponding  to  logical  operations  applied  to  them  (like  AND  and  OR)  in  a  very  compact  way. 
To  understand  the  real  practical  impact  of  this,  we  need  to  look  at  how  model  checking 
algorithms  work. 
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4.1.6  Model  checking  algorithms 

A  typical  model  checking  problem  consists  in  proving  that  a  given  property  is  satisfied  for  a 
certain  set  of  initial  states.  The  property  will  depend  on  the  type  of  logic  that  we  are  considering, 
but  typical  examples  are  whether  a  state  with  a  certain  characteristic  is  eventually  (or  always) 
reached,  or  that  a  certain  set  of  states  is  never  reached.  The  essence  of  model  checking 
algorithms  is  to  compute  the  set  of  states  that  satisfy  the  property  and  then  check  if  the  set  of 
initial  states  is  included  in  this  set.  The  advantage  of  using  BDDs  for  these  algorithms  is  that 
these  sets  are  never  defined  explicitly  (they  could  be  exponentially  large). 

As  an  example  of  the  basic  structure  of  these  algorithms,  let  us  consider  a  path  formula  in  CTL 
of  the  form  (O  U  'P),  which  means  that  for  all  paths  property  O  holds  until  property  ¥  becomes 
true.  The  goal  is  to  find  the  set  of  all  states  that  satisfy  this  formula.  The  algorithm  will  proceed 
iteratively  and  generate  a  nested  sequence  of  sets  that  satisfy  the  formula,  until  the  set  obtained 
does  not  change:  this  will  be  the  set  of  all  states  S  that  satisfy  the  formula,  and  then  we  can  easily 
check  if  the  set  of  initial  states  is  included  in  it. 

The  steps  in  such  an  iterative  algorithm  will  be  something  like  this: 

1 .  Let  Si  be  the  set  of  states  that  satisfy  'P  (sets  like  this  are  assumed  to  be  provided).  Then 
Si  ^  S,  since  all  states  that  satisfy  'P  trivially  satisfy  (O  U  'P). 

2.  Next,  find  the  set  of  states  that  both  satisfy  O,  and  such  that  there  is  a  transition  from  it  to 
a  state  in  Si.  Both  Si  and  the  set  of  states  that  satisfy  O  are  represented  by  BDDs.  The 
transition  relation  R,  being  a  subset  of  the  set  of  2N  bit  strings  can  also  be  represented  by 
a  BDD.  If  /i  is  the  Boolean  function  associated  with  Si,  and  x<t>  is  characteristic 
function  of  the  set  of  states  satisfying  O,  we  can  compute  the  Boolean  function  associated 
with  S2  ^  Si,  asf2ix)  =/i(x)  V  (/<j,(x)  A  3  x' .  (i?(x,  x')  A /i(x)).  The  key  point  is  that  all 
these  operations  can  be  computed  using  the  BDDs  associated  with  each  Boolean 
function. 

3.  Now  we  iterate  step  2,  looking  for  states  that  both  satisfy  O  and  have  a  transition  to  a 
state  in  S2. 

4.  Keep  iterating  until  the  BDD  associated  with  Sj  is  the  same  as  the  one  associated  with 
Sj+i. 

5.  Check  if  the  set  of  initial  states  is  included  in  Sj.  This  can  be  done  by  computing  the 
AND  of  the  BDD  associated  with  the  set  of  initial  states  and  the  one  associated  with  Sj 
and  checking  the  resulting  BDD  is  the  same  as  the  one  for  the  initial  states. 

As  we  can  see,  the  algorithm  looks  for  a  mathematical  object  (a  set)  that  is  extremal  in  the  sense 
that  is  the  larger  set  that  satisfies  the  property.  This  interpretation  gives  us  the  idea  that  the 
problem  may  be  recast  as  some  sort  of  optimization.  This  was  the  idea  put  forward  at  the 
Vanderbilt  meeting.  The  hope  was  that  such  an  optimization  problem  may  be  suitable  for 
implementation  with  DW2. 

4.1.7  Issues  with  casting  computation  of  extremal  BDD  as  optimization  problem 

In  order  to  cast  this  computation  of  an  extremal  BDD  as  a  combinatorial  optimization  problem 
suitable  for  implementation  with  DW2,  we  need  to  be  able  to  accomplish  two  tasks:  first,  encode 
BDDs  as  binary  strings  (the  space  over  which  DW2  performs  optimization),  and  second, 
construct  a  cost  function  that  whose  minimum  is  associated  with  the  extremal  BDD  we  are 
looking  for.  We  will  now  discuss  what  we  see  as  major  roadblocks  for  these  two  requirements. 
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BBD  encoding  as  binary  strings:  BDDs  are  used  to  eneode  Boolean  functions  on  N-bit  strings. 

Since  these  functions  are  uniquely  associated  the  subsets  of  N-bit  strings,  there  are  2  such 
functions.  However,  we  can  only  encode  a  space  of  size  2^  on  DW2.  So  here  is  the  first  major 
obstacle  to  optimizing  over  BDDs  using  DW2:  we  will  be  able  to  encode  all  BDDs  only  of 
systems  with  log  (Nqubits)  states,  and  since  the  number  of  qubits  is  currently  512  and  not 
expected  to  grow  exponentially,  this  approach  will  only  work  for  toy  models  and  will  not  scale 
for  even  moderately  sized  problems. 

A  way  around  this  issue  could  be  to  restrict  the  set  of  BDDs  over  which  the  optimization  is 
performed.  However,  there  is  no  clear  guidance  on  how  to  choose  this  subset.  Furthermore,  such 
a  subset  must  include  the  BDD  associated  with  the  solution  for  the  computation,  but  we  have  no 
way  of  knowing  if  that  will  hold.  If  we  proceed  anyway,  we  will  not  have  any  assurance  that  the 
BDD  obtained  has  any  of  the  properties  required  by  the  solution.  The  only  way  to  check  this  will 
be  to  perform  an  exhaustive  testing  of  such  a  BDD,  but  this  defeats  the  purpose  of  using  BDDs 
in  the  first  place. 

Cost  function  for  BDD  optimization:  even  if  we  assume  that  we  have  enough  qubits  to  encode  all 
possible  BDDs,  we  still  need  to  construct  a  suitable  cost  function  whose  minimum  is  associated 
with  the  extremal  BDD  sought.  From  the  algorithm  presented  above  it  is  clear  that  the  sequence 
of  BDD  generated  represent  subsets  of  the  state  space  with  increasing  size,  so  the  size  may  be  a 
candidate  for  a  cost  function.  But  it  is  also  clear  that  there  are  more  constraints  that  the  set  needs 
to  satisfy,  i.e.,  the  states  it  contains  need  to  satisfy  a  certain  formula  in  some  temporal  logic.  The 
parameters  at  our  disposal  when  constructing  a  cost  function  for  DW2  are  the  local  fields  and  the 
interaction  between  qubits.  This  results  in  a  quadratic  function,  and  it  is  not  clear  that  this  form 
has  the  power  to  encode  the  required  properties  of  the  states.  At  the  very  least,  this  mapping  (if 
possible)  would  require  a  lot  of  preprocessing  before  it  can  be  implemented  in  DW2,  and  this 
will  likely  erase  any  gains  produced  by  running  the  optimization  in  the  quantum  computer. 

4.1.8  Results  of  BDD  implementation  using  DW2 

As  discussed  above,  implementing  an  optimization  approach  to  solve  the  same  fix  point  problem 
that  is  usually  addressed  with  BDDs  has  several  obstacles.  It  may  be  the  case  that  these  obstacles 
can  be  overcome  in  some  cases  (for  example,  if  we  have  enough  knowledge  to  restrict  the  set  of 
BDDs  over  which  the  optimization  will  take  place),  but  at  that  point  in  the  research  we  believed 
that  the  cost  of  focusing  on  this  problem  and  identifying  favorable  instances  would  have  been  too 
high  for  the  expected  benefit.  Furthermore,  it  would  have  negatively  affected  the  research  on 
other  aspects  of  model  checking  that  seemed  more  likely  to  have  a  payoff.  It  was  then  decided 
not  to  pursue  this  avenue  of  research  as  part  of  the  QCHECK  project. 


4,2  Benchmarking  of  DW2  on  MAX-2-SAT  versus  MaxWalkSAT 

As  discussed  in  previous  quarterly  reports,  a  fair  comparison  of  the  performance  of  the  chip  is 
obtained  when  benchmarked  against  a  probabilistic  solver  like  MaxWalkSAT,  instead  of  an 
exact  solver  like  AK-MAXSAT.  An  exact  solver  will  typically  take  a  longer  time  to  run  since  it 
is  not  only  providing  the  solution,  but  also  a  guarantee  of  optimality.  On  the  other  hand,  the  D- 
Wave  processor  falls  under  the  class  of  probabilistic  solvers,  where  running  an  instance  will 
result  in  an  answer  that  corresponds  to  the  optimal  solution  with  probability  p<l.  Hence,  in  order 
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to  get  a  higher  confidence  in  the  result,  the  solver  needs  to  be  run  many  times,  and  this  must  be 
included  in  the  computation  of  the  runtime.  The  figure  of  merit  will  be  the  number  of  repetitions 
needed  to  find  the  best  answer  for  an  instance  with  a  probability  above  a  certain  threshold. 

To  run  this  benchmark  we  generated  random  instances  of  MAX-2-SAT  that  respected  the 
connectivity  of  the  Chimera  graph,  so  that  they  could  be  run  directly  on  the  D-Wave  processor. 
The  instances  had  N  variables  and  M=2N  clauses,  where  N  ranged  from  20  to  500  (in  steps  of 
20).  For  each  instance  we  used  an  exact  solver  to  determine  the  optimal  solution,  so  we  can  use  it 
to  compute  with  what  probability  the  D-Wave  processor  and  MaxWalkSAT  find  the  correct 
answer.  The  plot  below  shows  the  results  we  obtained  on  this  benchmark  using  a  2.5GHz 
desktop  processor  to  implement  MaxWalkSAT. 

4.2.1  Analysis  of  the  benchmarking  results 

The  results  of  Figure  7  show  the  comparison  between  the  performance  of  DW2  against  the 
MaxWalkSAT  algorithm  run  on  a  Mac  Pro  desktop  (2.6  GHz  processor).  We  plotted  the 
logarithm  of  the  runtime  (in  microseconds)  against  the  number  of  variables.  What  we  see  from 
this  experiment  is  that  the  DW2  processor  has  a  better  performance  for  problems  above  N=40 
variables.  This  is  an  encouraging  result,  although  it  should  be  fair  to  point  out  some  caveats. 
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Figure  7,  Benchmarking  results  for  MAX-2-SAT 

First,  the  MaxWalkSAT  algorithm  is  designed  to  tackle  general  MaxSAT  problems  and  no  steps 
were  taken  to  optimize  it  for  the  Chimera  connectivity  of  the  problems  in  the  ensemble  we 
tested.  In  some  sense,  one  can  see  this  issue  as  providing  some  kind  of  advantage  to  the  DW2 
processor,  since  the  problems  are  native  to  its  underlying  architecture.  This  is  an  issue  that 
should  always  be  kept  in  sight  when  interpreting  benchmarking  results  for  DW2. 
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Second,  we  have  only  run  these  instances  on  a  moderately  powerful  desktop  computer.  A  faster 
processor  will  certainly  increase  the  performance  of  MaxWalkSAT.  This  is  another  issue  that 
makes  comparing  the  DW2  quantum  device  with  a  classical  algorithm  somewhat  difficult.  We  do 
not  have  any  guidance  to  decide  against  which  classical  hardware  should  the  DW2  performance 
be  measured.  One  could  consider  the  cost  of  the  systems  (i.e.,  look  for  a  classical  computing 
system  of  similar  cost),  but  this  will  give  the  classical  approach  an  unfair  advantage,  given  that 
classical  hardware  and  algorithms  have  benefited  from  decades  of  research  and  investment. 
Quantum  systems  like  DW2  have  been  around  for  only  a  few  years  and  their  development  is  still 
at  the  very  early  stages. 

Another  issue  to  point  out  is  that  for  the  DW2  system  we  are  only  counting  the  time  required  to 
run  the  annealing  portion  of  the  optimization.  The  system  has  a  significant  time  overhead  that  is 
required  to  program  the  device  and  to  measure  the  results  at  the  end.  On  the  other  hand,  this 
overhead  is  constant  (and  being  improved  in  newer  designs)  and  will  play  less  of  a  factor  in 
future  generations  of  the  chip  having  a  larger  number  of  qubits. 

4.3  Implementation  of  a  heuristic  embedding  tool 

The  heuristic  embedding  tool  was  implemented  as  a  Matlab  program  that  communicated  with  the 
DW2  processor  through  a  special  function  call,  gathered  its  output  and  computed  the  gradient  of 
the  relative  entropy,  updated  the  Ising  model  and  submitted  it  again  to  the  DW2  processor.  The 
full  code  is  submitted  as  an  addendum  to  this  report.  We  will  now  discuss  the  main  components 
of  the  code. 

4.3.1  Code  structure 

The  main  code  is  the  Matlab  function  SEBREMf orQUBO  (SEBREM  stands  for  Sequential 
Embedding  By  Relative  Entropy  Minimization): 


o  Inputs:  1.  Qfull:  a  matrix  encoding  the  full  QUBO  problem  to  solve 

2.  beta:  a  parameter  that  characterizes  the  temperature  of  the  Gibbs 
distribution  associated  with  Qfull  (usually  taken  to  be  1). 

3.  Ni  terat  ions:  the  number  of  iterations  of  the  sequential  embedding. 

4.  EmbeddingFlag:  a  flag  that  determines  the  initial  embedding.  It  could 
take  three  values:  EmbeddingElag  =1,  applies  a  greedy  embedding  that  tries 
to  map  the  variables  to  the  qubits  in  a  way  that  preserves  the  couplers  with  the 
largest  absolute  value;  EmbeddingElag  =  2,  applies  a  modified  version  of 
the  greedy  embedding,  where  when  various  options  are  presented  when  choosing 
which  qubit  a  variable  will  be  assigned  to,  it  chooses  one  of  them  at  random  (the 
previous  embedding  would  choose  the  first  qubit  in  the  list);  EmbeddingElag 
=  3,  applies  a  randomized  direct  embedding,  i.e.,  it  randomly  assigns  variables 
to  qubits. 
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5.  Step:  controls  the  size  of  the  step  when  updating  the  {h,J)  parameters  in  the 
Ising  model. 

6.  Display:  (0  or  1)  turns  off  and  on  the  display  of  results  after  eaeh  iteration. 


o  Outputs:  1.  RelativeEntropy:  a  vector  of  values  of  the  relative  entropy  at  each 
iteration  step. 

2.  BestSolutions:a  matrix  whose  columns  are  the  best  solutions  found  at 
each  iteration  (a  solution  is  a  binary  string). 

3.  BestObj  ective:  a  vector  with  the  best  value  of  the  objective  function 
found  at  every  iteration. 

4.  BestSolutionFound:  a  binary  string  representing  the  best  solution  found 
over  all  iterations. 

5.  BestObj  ectiveFound:  the  value  of  the  objective  function  corresponding 
to  the  best  solution  found. 


The  eode  starts  by  fixing  the  values  of  certain  parameters  and  initializing  matrices  that  would 
store  the  results  generated  by  the  iterative  procedure.  Then  it  generates  the  initial  embedding 
(depending  on  the  value  of  EmbeddingFlag)  and  then  starts  the  iteration  that  construct  the 
sequential  embedding.  Eaeh  iteration  step  has  the  following  structure: 


1.  Solve  the  Ising  model  using  the  DW2  using  the  funetion  IsingConnectSolve,  whieh 
is  called  using  the  parameters  {h_chimera,J_chimera)  .  We  multiply  the  parameters  by 
0.5  in  order  to  inerease  the  number  of  samples  generated  (this  is  equivalent  to  raising  the 
temperature). 

2.  Extract  an  empirical  distribution  from  the  DW2  output,  i.e.,  the  solutions  sampled,  their 
frequency  and  their  assoeiated  energy. 

3.  Compute  the  objective  function  associated  with  the  matrix  Qfull  on  all  solutions 
sampled  in  the  previous  step. 

4.  Extract  the  solutions  that  minimize  the  value  of  the  objeetive  funetion. 

5.  Compute  the  relative  entropy. 

6.  Compute  the  gradient  of  the  relative  entropy. 

7.  Update  the  values  of  the  Ising  model  given  by  {h_chimera,J_chimera),  using  the  gradient 
information.  This  step  also  adds  another  heuristic  (implemented  inside  the  function 
UpdatelsingModel):  it  adds  to  the  Ising  model  another  Ising  model  that  has  the 
Chimera  connectivity  and  the  property  that  the  best  solution  seen  so  far  is  its  ground 
state.  The  idea  is  to  reward  direction  in  the  state  of  parameters  that  produced  good 
solutions.  This  heuristie  is  controlled  by  the  parameter  MIXFACTOR. 

8.  Information  is  displayed  (depending  on  the  value  of  the  Di  splay  flag)  and  saved. 
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9.  Go  back  to  step  1  until  Niterations  are  performed. 

4.3.2  Discussion  of  code  performanee 

The  main  objeetive  of  this  task  was  to  develop  a  tool  that  would  allow  us  to  solve  problems  with 
the  DW2  proeessor  that  would  not  fit  in  the  underlying  arehiteeture  given  by  the  Chimera  graph. 
In  that  capaeity,  the  code  generates  a  sequenee  of  Ising  models  that  respeet  the  Chimera  graph, 
with  parameters  that  are  adjusted  in  sueh  a  way  that  the  solutions  obtained  are  better  solutions  of 
the  original  non-Chimera  problem.  The  whole  method  ean  be  seen  as  an  optimization  in  the 
spaee  of  parameters  {h,J),  whieh  is  a  eontinuous,  bounded  subset  of  R"  (the  spaee  of  n-tuples  of 
real  numbers).  The  boundaries  are  given  by  the  limited  range  that  the  parameters  ean  take  when 
programming  the  deviee:  the  loeal  fields  lie  within  -2  and  2,  while  the  eouplers  take  values 
between  -1  and  1.  The  funetion  that  we  are  trying  to  optimize  in  this  spaee  is  the  relative  entropy, 
and  although  we  have  an  analytieal  expression  for  it,  the  function  is  not  convex  and  we  have  no 
guarantees  about  finding  the  absolute  minimum. 

It  is  then  hard  to  quantify  the  performanee  of  the  eode  in  general,  and  we  are  then  limited  to 
study  its  behavior  in  examples  of  interest.  The  basie  eharaeteristie  that  we  would  want  the  eode 
to  have,  is  that  it  would  improve  the  solutions  provided  by  the  initial  embedding.  The  initial 
embedding  is  in  some  sense  the  best  we  ean  do  to  solve,  in  one  shot,  a  problem  that  does  not  fit 
the  proeessor’ s  eonneetivity.  Different  strategies  for  this  embedding  will  result  in  different 
quality  of  solutions.  We  eonsidered  three  types  of  initial  embeddings:  two  that  tried  to  eapture 
some  of  the  strueture  of  the  problem  (greedy  embeddings),  and  one  that  essentially  mapped  the 
problem  into  the  proeessor  in  a  random  way.  An  interesting  feature  we  found  is  that,  even  though 
the  greedy  strategies  provided  better  initial  solutions  than  the  randomized  embedding,  by 
implementing  the  sequential  embedding  we  were  always  able  to  generate  better  solutions  after  a 
few  iterations,  irrespeetive  of  the  initial  embedding  we  ehose.  Given  that  the  greedy  embeddings 
required  some  non  trivial  preproeessing  (i.e.,  finding  sueh  embedding)  that  required  extra 
eomputational  resources,  we  realized  that  a  randomized  embedding  would  provide  good 
solutions  without  the  upfront  eomputational  eost. 

The  time  performanee  of  the  sequential  embedding  depends  on  the  size  of  the  problem  being 
eonsidered,  mainly  beeause  the  eomputation  of  the  gradient  ean  be  expensive  if  eare  is  not  taken 
to  eode  it  effieiently.  We  implemented  the  eomputation  in  Matlab,  and  we  took  eare  of 
veetorizing  the  ealeulations  as  mueh  as  possible  in  order  to  avoid  any  loops  (that  are  notoriously 
slow  in  Matlab).  For  problems  of  around  500  variables,  the  gradient  computation  step  would  take 
only  a  few  seeonds.  We  found  that  in  most  problems,  only  a  few  tens  of  iterations  were  enough 
to  reaeh  a  point  where  the  solutions  would  no  longer  improve  (we  suspeet  we  were  either 
reaehing  a  loeal  minima  of  the  relative  entropy,  or  that  the  approximations  we  made  were  no 
longer  valid).  So  in  summary,  for  a  500  variable  problem,  a  run  time  of  a  few  minutes  will  be  the 
most  we  needed  to  run  in  order  to  find  the  best  solutions  this  method  would  provide.  Applying 
this  method  to  future  generations  of  the  proeessor,  that  would  support  more  than  a  thousand 
variables,  would  eertainly  benefit  from  speeding  up  the  gradient  computation  by  using  faster 
languages. 

The  algorithm  has  several  parameters  that  can  be  used  to  try  to  improve  convergence  and 
behavior.  In  this  projeet  we  did  not  have  the  time  required  to  analyze  them  in  more  detail  and  try 
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to  find  their  optimal  values.  The  values  we  used  were  bom  out  of  trial  and  error,  where  the 
driving  feature  was  that  the  code  generated  improving  solutions  in  a  reasonable  amount  of  time. 

In  summary,  the  tool  we  developed  succeeded  in  allowing  us  to  produce  good  solutions  for 
optimization  problems  that  did  not  fit  into  the  processor’s  connectivity.  These  solutions  were 
better  than  the  ones  obtained  using  single-shot  approximate  embeddings.  We  can  always  take 
any  approximate  embedding  as  the  initial  embedding,  and  we  believe  that  the  algorithm 
implemented  here  will  always  improve  the  solutions  initially  obtained. 

4,4  Integration  of  DW2  into  CEGAR  loop 

As  discussed  in  previous  sections,  the  central  goal  of  this  project  was  to  integrate  the  capabilities 
of  the  DW2  processor  into  the  CEGAR  framework  of  model  checking.  We  identified  two  key 
steps  in  the  CEGAR  loop  where  combinatorial  optimization  problems  needed  to  be  solved  in 
order  to  proceed; 

i.  Solving  a  SAT  problem  to  verify  the  validity  of  an  abstract  counterexample 

ii.  Solving  an  lEP  problem  to  find  the  smallest  increase  of  the  abstraction  in  order  to  get  rid 
of  a  spurious  abstract  counterexample 


The  first  of  the  two  tasks,  at  least  in  the  form  we  were  able  to  cast  it,  turned  out  not  to  be  a  good 
fit  for  the  capabilities  of  the  DW2  processor.  The  second  task  provided  much  more  promising 
results. 

4.4.1  Verifying  the  validity  of  abstract  counterexamples 

Given  an  abstract  counterexample,  checking  its  validity  consists  in  verifying  if  the  corresponding 
trace  in  the  original  system  also  provides  a  counter  example  to  the  property  we  are  trying  to 
prove.  This  problem  can  be  formally  reduced  to  checking  the  satisfiability  of  a  Boolean  formula 
that  encodes  the  existence  (or  non  existence)  of  the  required  trace  in  the  original  system.  In  order 
to  study  this  task,  we  generated  abstract  problems  in  New  Symbolic  Model  Checker  for  System 
Verification  (NuSMV),  we  checked  the  abstract  models,  and  reformulated  the  counterexample 
checking  problem  in  Einear  Temporal  Logic  (LTL).  We  wrote  a  program  that  read  abstract 
counter-examples,  and  generated  from  them  LTL  formulas  that  would  be  satisfiable  if  and  only  if 
the  abstract  counterexamples  were  valid,  i.e.,  they  were  associated  with  actual  counterexamples 
in  the  original  system.  Then  we  exported  the  resulting  SAT  problem  using  a  canonical  format, 
converted  them  into  QUBO  problems  and  tested  on  DW2. 

Our  first  experiments  were  with  the  problem  of  testing  abstract  counter-examples  for  soundness. 
We  then  translated  these  LTL  formulas,  and  the  system  models,  into  satisfiability  problems, 
using  a  standard  format  known  as  DIMACS.  These  SAT  problems  were  then  be  reformulated 
and  submitted  to  the  DW2  to  be  solved. 

We  worked  with  the  NuSMV  system  developed  by  Londazione  Bruno  Kessler  (LBK),  CMU,  and 
the  Universities  of  Trento  and  Genova  (http://nusmv.fbk.eu).  NuSMV  is  the  next  generation, 
open  source  version  of  the  original  Symbolic  Model  Checker  for  System  Verification  (SMV) 
model-checking  tool.  SMV  was  the  first  BDD-based  model-checker,  developed  by  McMillan  at 
CMU.  NuSMV  extends  the  original  SMV,  and  offers  both  BDD-based  model-checking  and 
SAT-based  bounded  model-checking.  NuSMV  provided  the  model-checking  of  the  abstract 
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models,  and  translated  our  LTL  elaims  into  SAT  problems  (NuSMV  provides  SAT  problem 
translation  to  support  external  SAT  tools  for  BMC). 

For  our  experiments  with  soundness  eheeking,  we  manually  abstraeted  various  models  from  the 
NuSMV  distribution.  These  models  included  some  properties  to  check,  and  we  added  other 
properties  of  our  own.  We  added  some  of  our  own  properties  since  most  of  the  properties  in  the 
distribution  were  properties  that  the  models  satisfied,  so  we  needed  to  have  additional  cases 
covering  unsafe  models. 

To  generate  the  abstract  counter-examples,  we  used  NuSMV’s  BDD-based  model-checking  on 
the  manually  abstracted  models.  NuSMV  can  be  configured  to  write  counter-examples  in  an 
Extended  Markup  Language  (XML)  form  that  is  very  easy  to  parse.  We  are  grateful  to  LBK 
personnel  for  helping  us  understand  the  XML  format,  and  for  responding  to  our  issues  with  the 
format. 

Reformulation:  The  key  to  our  checking  process  is  reformulating  the  validity  checking  problem 
as  an  LTL  claim,  since  this  enables  automatic  generation  of  the  SAT  problem.  The  abstract 
counter-example  is  a  sequence  of  abstract  states,  asi,  each  of  which  is  an  assignment  of  values  to 
propositional  variables.  Since  these  are  abstract  states,  the  assignments  are  to  only  a  subset  of 
the  propositional  variables  of  the  full  model.  Again,  since  this  is  an  abstraction,  a  single  abstract 
transition  may  correspond  to  multiple  concrete  transitions.  So  for  an  abstract  counter-example 
aso,asi...aSn  to  be  valid,  there  must  be  a  valid  sequence  of  concrete  states  so...Si...Sm  such  that  aso 
^  So,  as  I  ^  Si,, „  asn^  Sm-  Note  that  the  subset  relation  is  equivalent  to  saying  that  the  concrete 
state  entails  the  abstract  state.  We  may  formulate  this  in  LTL  as  follows: 

asoX(asi  X(as2  ...  X as„)))  (Property  1)  (24) 

That  is,  we  begin  in  a  state  satisfying  aso,  eventually  we  reach  a  state  satisfying  asj,  then 
eventually  we  reach  a  state  satisfying  as2,  and  so  on,  until  we  reach  a  state  satisfying  as-a,  all  in 
the  context  of  the  full,  original,  concrete  model. 

The  above  query  cannot  be  used  as  written  to  check  validity,  however.  As  an  LTL  property,  it  is 
effectively  claiming  that  in  all  runs  of  the  model,  we  start  in  a  state  satisfying  aso,  reach  a  state 
satisfying  asi,  etc.  To  use  this  formulation,  we  must  invert  it,  and  charge  the  solver  to  prove  to 
us  that  it  is  impossible  to  satisfy  Property  1 .  If  the  negation  of  Property  1  is  valid,  then  the 
abstract  counter-example  is  invalid: 

—asoX(as ;  X(as2  ...  X as„)))  (Property  2)  (25) 

Implementation:  To  support  these  experiments,  we  wrote  a  program  that  parsed  the  XML 
formatted  counterexamples,  and  generated  from  them  LTL  claims  in  the  form  of  Property  2, 
expressed  in  NuSMV’s  input  language.  These  LTL  claims,  together  with  the  original  (concrete) 
system  model,  were  then  submitted  to  NuSMV,  and  bounded  model  checking  SAT  models  in 
DIMACS  notation  were  extracted.  These  DIMACS  problems  were  then  translated  into  problems 
for  the  D-Wave  system. 

Unfortunately,  experience  has  shown  that  the  D-Wave  quantum  adiabatic  optimization  process  is 
not  a  good  choice  for  this  kind  of  decision  problem.  The  D-wave  is  well  suited  to  approximate 
optimization  problems,  where  best  effort  is  what  matters.  But  in  the  case  of  decision  problems,  a 
best  effort  that  misses  is  not  a  useful  approximation;  it  is  simply  wrong.  The  CEGAR  algorithm 
shows  that  approximate  solutions  that  can  be  wrong  are  sometimes  acceptable,  if  the 
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approximations  are  safe  (conservative),  and  if  there  is  a  means  to  refine  the  computation.  So  it  is 
possible  that  a  revision  of  our  approach  here,  where  the  D-Wave  would  only  give  erroneous 
results  that  mistakenly  concluded  that  the  counter-example  was  valid,  and  never  give  a  false 
invalid  result,  would  be  useful.  Such  a  technique  would  never  cause  us  to  mistakenly  certify  a 
system  as  safe  when  it  wasn’t,  and  would  never  send  us  on  a  wild-goose  chase  to  refine  an 
already  sound  counter-example. 

4.4.2  Refining  the  model 

The  second  part  of  the  QCHECK  project  involved  using  the  D-wave  to  solve  the  optimization 
problems  arising  in  the  refinement  part  of  the  CEGAR  algorithm.  Recall  that  this  involves 
finding  an  optimal  refinement  of  the  model.  In  the  specific  version  of  CEGAR  that  we  used,  this 
was,  specifically,  finding  a  minimal  set  of  propositional  variables,  previously  abstracted  away,  to 
add  to  the  abstract  model.  This  is  formulated  as  an  ILP,  choosing  a  set  of  variables  that  are 
sufficient  to  create  a  new  abstraction  in  which  the  dead  end  states  and  the  bad  states  of  an  invalid 
counter-example,  are  guaranteed  to  be  in  different  abstract  states. 

Our  work  in  this  area  was  more  successful  than  our  work  on  checking  abstract  counter-examples. 
This  was  primarily  due  to  the  fact  that  the  abstraction-refinement  is  a  better  fit  to  the  D-wave ’s 
capabilities.  In  particular,  this  problem  is  an  optimization  problem,  rather  than  a  decision 
problem,  and  if  the  solution  is  less  than  perfectly  optimal,  it  simply  costs  us  more  work  in  the 
model-checking  phases  of  the  CEGAR  process:  a  subop timal  answer  can  never  lead  to  an 
unsound  conclusion. 

Our  original  plan  was  to  find  an  off-the-shelf  model-checker  with  an  implementation  of  the 
CEGAR  algorithm,  from  which  we  could  extract  ILPs  or,  as  a  second  choice,  a  model  checker  in 
which  we  could  perform  steps  of  the  algorithm  and  from  which  we  could  extract  partial  results. 
Unfortunately,  our  attempts  to  find  such  a  tool  were  not  successful.  We  contacted  the  authors  of 
the  original  CEGAR  paper,  and  while  they  offered  a  number  of  very  helpful  suggestions,  the 
code  for  the  original  implementation  had  been  lost.  We  consulted  EBK  about  NuSMV,  and 
although  its  successor,  NuXmv,  will  contain  an  implementation  of  the  CEGAR  algorithm,  it  was 
not  in  a  condition  for  release  to  us.  ^  We  also  investigated  whether  it  would  be  possible  to  work 
the  CEGAR  algorithm  “around”  a  tool,  in  the  way  we  were  able  to  do  our  counter-example 
checking  “around”  the  NuSMV  tool.  We  concluded  that  this  would  not  be  possible,  since  the 
data  structures  needed  to  identify  the  dead  end  and  bad  states  were  not  visible  to  the  user  of 
NuSMV.  We  also  investigated  whether  it  would  be  possible  to  implement  the  CEGAR 
algorithm  by  scripting  NuSMV.  This  technique  seemed  very  promising:  the  NuSMV  code  is 
well  structured,  and  we  were  able  to  construct  NuSMV  data  structures  and  exercise  some  of  its 
functions  through  a  foreign  function  interface  from  the  Common  Eisp  programming  language. 
Unfortunately,  we  concluded  that  for  this  small  program,  we  didn’t  have  the  resources  to  dig 
deeply  into  the  NuSMV  sources.  Scripting  NuSMV  would  be  a  promising  direction  to  take  given 
a  larger-scale  program. 


^  It  has  since  been  released,  but  so  far  only  in  closed-source  form,  so  would  not  be  suitable  for 
our  purposes. 
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After  we  had  to  abandon  the  plan  of  using  an  off-the-shelf  tool,  our  next  step  was  to  manually  go 
through  the  steps  of  the  CEGAR  algorithm,  using  two  BDD  “desk  caleulators,”  to  eheck  the 
abstraetions,  and  to  identify  dead  ends  and  bad  states.  In  particular,  we  used  both  the  iben  and 
bddc^  tools.  We  learned  a  number  of  valuable  lessons  about  the  structure  of  models  that  would 
give  rise  to  valid  and  invalid  abstract  counter-examples,  and  designed  some  small  problems. 
After  developing  a  small  number  of  small  problems,  we  were  unable  to  make  larger  problems: 
the  need  to  work  directly  with  propositional  logic  was  too  burdensome  and  inefficient. 

At  this  point  we  were  nearing  the  end  of  the  program,  and  realized  that  we  needed  our  own 
implementation  of  the  CEGAR  algorithm,  so  that  we  could  run  examples  more  easily.  We  were 
fortunate  to  find  that  the  CU  Decision  Diagram  (CUDD)  open  source  Ordered  Reduced  Binary 
Decision  Diagram  (OBDD)  library,"^  developed  by  Eabio  Somenzi  at  Colorado  University,  offers 
a  Perl  application  programming  interface  (API).  Using  this  API  we  could  rapidly  build  our  own 
implementation  around  the  BDD  operations,  developing  interactively. 

The  CUDD  library  provided  one  of  the  pieces  of  the  solution:  it  remained  to  make  or  choose  an 
input  language.  We  considered  the  SMV  language,  but  concluded  that  writing  a  parser  would 
require  too  much  effort.  We  had  consulted  Ofer  Strichman,  one  of  the  authors  of  the  original 
CEGAR  paper,  and  he  suggested  that  we  use  the  AIGER  -  And  Inverter  Graph  -  notation.^  The 
AIGER  format  offered  a  number  of  advantages:  (I)  it  is  very  easy  to  parse;  (2)  it  is  used  for 
hardware  model-checking  competitions,  so  there  are  a  large  number  of  preexisting  models 
available;  (3)  it  maps  nicely  to  the  abstraction  framework  of  CEGAR.  With  respect  to  point  (3), 
in  addition  to  AND  gates  and  Inverters,  AIGER  models  contain  latches.  A  CEGAR-compatible 
abstraction  technique  is  to  abstract  away  latches  by  treating  them  as  if  they  are  inputs.  Einally, 
point  (4),  there  are  existing  BMC  tools  for  AIGER  models  that  we  could  (and  did)  use  to  check 
our  results  when  building  and  debugging  our  CEGAR  implementation. 

In  practice,  the  AIGER  notation  did  not  provide  a  perfect  input  solution:  although  it  is  higher 
level  than  the  propositional  logic  that  we  used  when  working  with  the  BDD  calculators,  it  was 
still  too  low  level.  In  particular,  AIGER  offered  no  way  to  capture  higher  level  building  blocks 
such  as,  for  example  OR-gates,  XOR-gates,  half-adders,  etc.  In  order  to  build  models  that  would 
offer  interesting  lEPs,  we  needed  to  be  able  to  reuse  model  components.  To  support  such  reuse, 
we  added  a  simple  macro-language,  laig  (for  Eisp-fiavored  AIG),  and  wrote  a  simple  macro¬ 
expansion  facility.  We  discuss  this  in  more  detail  below. 

Three  other  issues  remained.  Eirst,  although  the  laig  notation  made  it  possible  to  reuse  model 
components,  allowing  us  to  make  bigger  and  more  interesting  circuits,  it  did  nothing  to  help 
formulate  the  properties  we  needed  to  check.  In  our  successful  test  cases,  we  were  checking 
properties  of  the  form  “in  a  valid  trace  it  is  impossible  to  reach  a  state  satisfying  P.  ”  The  proviso 
“in  a  valid  trace”  was  necessary  in  order  to  capture  bounding  assumptions  such  as  “there  are  no 


2 

http://sourceforge.net/projects/iben/ 

http://www-verimag.imag.fr/PEOPEE/Pascal.Raymond/tools/bddc-manual/bddc-manual- 

pages.html 

^  http://vlsi.colorado.edu/~fabio/CUDD/ 

^  http://fmv.jku.at/aiger/ 
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more  than  two  sensor  failures”  in  our  example  avionies  model.  In  practiee,  assumptions  like 
“valid  trace”  unpacked  into  temporal  properties  which  were  difficult  to  formulate  in  AIGER  and 
which  laig  notation  did  not  support  as  well  as  possible.  To  capture  temporal  properties  in 
AIGER  notation  involved  adding  latches  that  captured  property  components  so  that  temporal 
assertions  could  be  evaluated.  While  laig  allowed  us  to  capture  repeated  components  of  such 
translations,  an  automatic  translation  of  temporal  properties  would  have  been  helpful.  Another 
problem  was  that  the  existing  AIGER  models  were  not  very  well  documented,  so  they  didn’t 
provide  all  the  CEGAR  examples  we  had  wanted.  Nevertheless,  the  availability  of  examples, 
with  gold  standard  computation  results  (from  the  BMC  checker),  was  immensely  useful  in 
development  and  debugging.  Einally,  the  abstraction  relationship  of  treating  latches  as  inputs 
turned  out  to  be  more  subtle  than  we  had  originally  anticipated;  different  treatments  of 
abstraction  in  model-checking  used  subtly  different  definitions,  so  that  getting  consistent  results 
required  substantial  debugging  and  rework. 

4.4.3  CEGAR  Implementation 

To  summarize  the  implementation,  it  contained  the  following  components: 

Input: 

•  The  Eisp-style  AIGER  language  (laig)  and  translator.  Input  models  are  formulated 
as  laig  files,  and  translated  into 

•  Single-output  AIGER  models.  Single-output  models  are  used  to  capture  reachability 
checks  in  AIGER  models.  A  model  formulated  in  this  way  is  safe  if  it  can  never 
output  1. 

•  A  single-output  AIGER  model  is  read  by  our  AIGER  parser,  and  translated  into 
internal  data  structures. 

CEGAR  Process: 

•  An  initial  abstraction  is  computed,  which  hides  all  of  the  latches  that  are  not  direct 
components  of  the  single  AIGER  output. 

•  EOOP 

o  Our  solver  checks  reachability  of  the  unsafe  state  in  the  abstract  model.  This  is 
done  using  a  search,  implemented  as  a  Perl  loop  (based  upon  code  from  Ed 
Clarke’s  model-checking  text),  whose  primitive  computations  are  done  using  the 
CUDD  BDD  library. 

o  If  the  unsafe  state  is  not  reachable,  we  return  “safe”  and  exit. 

o  If  the  unsafe  state  is  reachable,  we  compute  an  abstract  counterexample. 

o  We  check  the  abstract  counterexample  for  validity.  This  is  done  according  to  the 
CEGAR  algorithm,  using  a  loop  over  the  images  computed  in  the  reachability 
search. 

o  If  the  abstract  counterexample  is  valid,  we  return  “unsafe”  and  exit. 
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o  If  the  abstract  counterexample  is  invalid,  we  identify  the  dead  end  and  bad  states. 
We  compute  an  ILP,  each  of  whose  rows  captures  the  difference  between  one 
dead  state,  d,  and  one  bad  state,  b.  The  columns  of  the  ILP  are  latches  that  are 
hidden,  and  each  row  assigns  1  to  propositional  variables  whose  values  differ  in  b 
and  d,  and  zero  to  all  other  entries. 

o  The  CEGAR  algorithm  publishes  the  ILP  to  the  D-wave  solver,  and  reads  a 
solution  to  the  ILP.  The  solution  to  the  ILP  specifies  a  set  of  latches  that  should 
no  longer  be  hidden. 

o  We  compute  a  new  abstract  model,  by  “unhiding”  the  set  of  latches  indicated  by 
the  D-Wave  solver  and  go  to  LOOP. 


We  were  pleased  to  find  that  it  was  very  easy  to  interface  the  CEGAR  loop  with  the  existing  D- 
Wave  solver.  The  output  format  was  very  easy  to  print  and  to  parse,  as  was  the  format  of  the 
answers  (a  list  of  latch  indices).  This  holds  promise  for  further  work  on  the  CEGAR  algorithm, 
and  for  other  applications  in  which  an  outer  loop  is  to  be  wrapped  around  the  D-Wave  solver. 


4.4.4  lEP  Problems 


Particular  structures  in  model-checking  problems  pose  interesting  challenges  for  abstraction 
refinement.  Or,  put  differently,  many  spurious  counter-examples  can  easily  be  eliminated  by 
adding  one  or  two  new  latches.  More  challenging  problems  arise  in  the  presence  of  two  features. 
Eirst,  there  must  be  multiple  different  concrete  paths  that  correspond  to  a  single  spurious  abstract 
counterexample.  This  arises  when  there  are  functions  in  the  counterexample  that  involve  a  wide 
variety  of  variables  in  a  non-trivial  Boolean  function.  The  second  condition  is  that  these  Boolean 
functions  must  not  include  intermediate  values  that  are  computed  and  latched.  If  there  is  such 
structure,  then  the  CEGAR  algorithm  will  simply  identify  whether  or  not  the  set  of  values  so 
latched  are  consistent  and  either  confirm  the  abstract  counter-example  or  reject  it,  without  need 
for  complex  reasoning  (this  is  the  strength  of  the  CEGAR  algorithm).  We  conjecture  that  the 
challenging  cases  for  CEGAR  refinement  are  characteristic  of  models  of  complex  circuits.  They 
may  be  less  likely  to  arise  in  models  of  software,  where  storage  of  intermediate  results  is  more 
common.  As  the  D-Wave  system  gets  more  powerful,  and  capable  of  handling  larger  and  larger 
models,  we  believe  that  we  will  see  more  non-trivial  ILPs  in  CEGAR  refinement. 

4,5  Implementation  of  CEGAR  based  model  checking  example 

We  describe  here  an  example  run  of  the  CEGAR  algorithm  with  the  D-Wave  solver  providing 
answers  to  the  ILPs.^  The  test  model  ^  was  inspired  by  an  avionics  example  problem  (the 


^  The  CEGAR  program  can  be  run  without  the  D-Wave  in  the  loop,  simply  picking  a  single 
implicated  variable  to  “unhide”  at  each  iteration.  Of  course,  this  can  perform  arbitrarily  badly;  it 
is  only  of  use  on  small  test  runs,  to  debug  other  parts  of  the  CEGAR  loop. 
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Rockwell  Flight  Control  System  (FCS)  5000  problem  [17]),  but  does  not  eorrespond  to  any  real 
hardware  or  software.  The  property  eheeked  is  a  mutual  exelusion  property:  there  are  two 
eontrol  modes,  “left”  and  “right,”  and  the  invariant  is  that  the  system  should  never  have  both 
eontrol  modes  aetive  simultaneously. 

4.5.1  Verihe  ation  summary 

Initially,  the  model  is  almost  eompletely  abstraet.  Reeall  that  the  property  to  be  eheeked  is  of  the 
form  “in  a  valid  traee,  the  left  and  right  eontrol  modes  eannot  both  be  engaged.”  The  only  lateh 
involved  in  this  claim  is  the  “traee  valid  lateh.”  This  allows  for  a  trivial  eounterexample,  since 
“left  eontrol  mode  on”  and  “right  eontrol  mode  on,”  are  essentially  treated  as  inputs.  This  is  an 
invalid  eounterexample,  for  reasons  having  to  do  with  system  funetion  —  there  is  substantial 
logie  behind  whether  or  not  the  eontrol  modes  are  engaged  -  and  the  logic  of  the  property  — 
there  is  substantial  logic  in  the  memory  of  traeking  whether  a  given  traee  is  valid  or  not.  The 
first  several  refinements  do  not  require  any  optimization:  they  are  obvious  refinements  that 
gradually  add  logie  for  the  full  temporal  ehain  needed  to  deeide  traee  validity. 

The  first  interesting  ILP  arises  when  the  model  has  been  elaborated  around  the  temporal  logic  of 
deciding  trace  validity.  There  is  now  a  spurious  eounterexample  to  the  effeet  of  “turn  the  right 
eontroher  on,  turn  the  left  eontroher  on,  wait  for  three  eyeles,  then  turn  the  valid  lateh  on.”  The 
model  actually  enforces  that  the  control  modes  ean  only  be  engaged  if  sensors  detect  the  right 
eondition  and  enough  of  the  redundant  sensors  vote  to  enable  it.  This  meets  the  eriterion  above 
for  an  interesting  ILP:  there  are  many  eombinations  of  the  sensor  latehes  that  eould  lead  to  a 
sueeessful  vote  to  engage  the  sensor  modes,  and  the  ILP  solver  must  ehoose  a  minimal  eovering 
set. 

The  seeond,  and  last,  interesting  ILPs  arise  from  a  spurious  counterexample  in  whieh  the  sensor 
inputs  are  arranged  in  sueh  a  way  that  both  the  left  and  right  eontrol  modes  ean  be 
simultaneously  engaged.  This  is  a  spurious  eounterexample  beeause  of  the  traee  validity 
eheeking:  it  turns  out  that  such  a  constellation  of  sensor  inputs  eannot  happen  in  valid  traees. 
Note,  however,  that  we  don’t  simply  eliminate  all  ineonsistent  sensor  value  eombinations:  we 
simply  stipulate  a  limit  on  the  number  of  sensors  that  ean  simultaneously  give  the  wrong 
deteetions.  The  traee  validity  eonstraints  eapture  the  notion  that,  e.g.,  the  aireraft  eannot  be 
taking  off  and  landing  at  the  same  time.  In  solving  the  final  interesting  ILP,  the  solver  identifies 
a  set  of  variables  that  is  suffieient  to  eapture  enough  of  the  validity- eheeking  logic  to  eliminate 
these  ineonsistent  sensor  eombinations. 

This  summary  has  skipped  over  the  refinement  steps  that  do  not  require  the  intervention  of  the 
ILP  solver.  There  are  many  sueh  steps,  where  identifying  a  set  of  latehes  to  add  to  the  model  is 
trivial.  This  arises  when,  e.g.,  the  same  lateh  appears  in  ah  rows  of  the  ILP.  In  those  eases,  the 
CEGAR  loop  does  not  invoke  the  D-Wave.  We  give  the  full  transeript  below. 

4.5.2  Detailed  transeript 

Initially,  the  model  is  almost  eompletely  abstract.  The  only  un-hidden  variable  is  the  one 
eorresponding  to  “the  traee  is  valid”  (literal  52).  Unsurprisingly,  this  leads  to  a  trivial 


’  We  give  the  laig  souree  as  an  appendix. 
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counterexample:  we  just  turn  on  the  right  and  left  eontrollers  (sinee  they  are  uneonstrained  in  this 
abstraet  model). 

This  abstraetion  is  too  simple,  for  two  different  reasons:  The  first  is  the  interesting  reason,  whieh 
is  that  the  eounterexample  must  take  into  aeeount  the  surrounding  environment,  as  redundantly 
sensed.  The  less  interesting  reason  is  that  there  is  a  ehain  of  latehes  for  the  on/off  state  of  the  left 
and  right  eontrollers;  this  is  neeessary  to  synehronize  the  state  of  the  system  with  the  state  of  the 
part  of  the  eircuit  that  traeks  validity  of  the  traee.  A  traee  that’s  invalid  (e.g.,  we  are  both  at  and 
not  at  an  airport,  simultaneously)  must  be  weeded  out. 

To  be  preeise,  we  are  verifying  that  we  never  enter  the  bad  state  of  “leftS  is  on,”  “rights  is  on,” 
and  “the  traee  is  valid.” 

The  first  several  steps  of  the  CEGAR  refinement  proeess  introduee,  one  by  one,  the  latehes  in  the 
three-step  delay  chains  from  “rightl”  to  “rightS”  and  “leftl”  to  “leftS.” 

Made  new  Aiger  model  with  9  inputs,  40  latches,  1  outputs,  and  42 
and  gates. 

Reading  40  latches. 

Reading  1  outputs. 

Reading  42  and  gates. 

Done  reading  file  and  building  model. 

Checking  abstract  model,  hidden  literals  are:  20  22  24  26  28  30  32 

34  36  38  40  42  44  46  48  50  54  56  58  60  62  64  66  68  70  72  74  76  78  80 

82  84  86  88  90  92  94 

Making  transition  relation 

Done  making  transition  relation 

Preparing  for  model-checking. 

Starting  fixed  point  computation. 

Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Goal  reached 
Goal  Reached  in  2  steps 
Counterexample : 

0 

Inputs:  111111111 

Latches:  0000000000000000000000000000000000000000 
Output :  0 

1 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIOIIIIIIIIIIIIIIIIIIIOO 
Output :  0 

2 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII 
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Output :  1 

Abstract  counterexample: 

000 

100 

111 

Abstract  counterexample  is  invalid  at  step  2 . 

ILP  is: 

92  94 

Splitting  on  literal  92 

The  first  refinement  is  to  add  one  of  the  two  latches  (literals  92  &  94),  representing  the  state  of 
the  left  and  right  controllers.  In  this  case,  the  algorithm  adds  the  literal  92  (“leftS”).  This  does  not 
require  the  intervention  of  the  D-Wave  solver;  the  ILP  is  trivial  (the  format  “92  94”  represents  a 
linear  constraint  in  which  the  sum  of  literals  92  and  94  must  be  greater  or  equal  than  1). 

Checking  for  goal 
Goal  reached 
Goal  Reached  in  3  steps 
Counterexample : 

0 

Inputs:  111111111 

Latches:  0000000000000000000000000000000000000000 
Output :  0 

1 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIOIIIIIIIIIIIIIIIIIOIOO 
Output :  0 

2 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIOIIIIIIIIIIIIIIIIIIIOI 
Output :  0 
3 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII 
Output :  I 

Abstract  counterexample: 

0000 

1000 

1101 

1111 

Abstract  counterexample  is  invalid  at  step  2 . 

ILP  is: 

88  94 

Splitting  on  literal  88 


The  next  counterexample  is  almost  as  trivial.  The  solver  realizes  that  the  left  channel  cannot 
simply  be  turned  on  instantaneously,  but  still  generates  a  very  quick  path  to  failure  (only  three 
steps  long),  taking  into  account  a  delay  in  turning  the  left  channel  on.  Note  that  this  delay  comes 
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partly  because  there  is  a  delay  from  the  logic  of  checking  the  switch  and  the  state-sensing  logic, 
but  also  because  the  model  has  a  three-step  delay  to  compute  the  validity  of  the  trace. 

This  counterexample  is  shown  to  be  spurious:  the  abstract  model  now  takes  into  account  the 
delay  in  turning  on  the  left  controller,  but  not  yet  the  right  controller.  The  culprit  for  this  three- 
step  counterexample  is  that  there  must  be  a  progression  (at  least)  through  latches  “right2”  (88) 
and  “rights”  (94).  The  solver  has  another  trivial  refinement  problem,  to  choose  between  88  and 
94,  which  it  does  by  choosing  in  numerical  order  (again,  no  need  for  the  ILP  solver;  there  is  only 
one  row  to  the  ILP). 


Starting  fixed  point  computation. 

Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Goal  reached 
Goal  Reached  in  4  steps 
Counterexample : 

0 

Inputs:  111111111 

Latches:  0000000000000000000000000000000000000000 
Output :  0 

1 

Inputs:  111111111 

Latches :  1111111111111111110111111111111111010100 
Output :  0 

2 

Inputs:  111111111 

Latches :  1111111111111111110111111111111111110101 
Output :  0 

3 

Inputs:  111111111 

Latches :  1111111111111111110111111111111111111101 
Output :  0 

4 

Inputs:  111111111 

Latches :  1111111111111111111111111111111111111111 
Output :  1 

Abstract  counterexample: 

00000 
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10000 

11001 

11101 

mil 

Abstract  counterexample  is  invalid  at  step  2 . 

ILP  is: 

82  94 

Splitting  on  literal  82 

Checking  abstract  model,  hidden  literals  are:  20  22  24  26  28  30  32 
34  36  38  40  42  44  46  48  50  54  56  58  60  62  64  66  68  70  72  74  76  78  80 
84  86  90  94 

Making  transition  relation 
Done  making  transition  relation 
Preparing  for  model-checking. 

Starting  fixed  point  computation. 

Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Goal  reached 
Goal  Reached  in  5  steps 
Counterexample : 

0 

Inputs:  111111111 

Latches:  0000000000000000000000000000000000000000 
Output :  0 

1 

Inputs:  111111111 

Latches :  1111111111111111110111111111111011010100 
Output :  0 

2 

Inputs:  111111111 

Latches :  1111111111111111110111111111111111010101 
Output :  0 
3 

Inputs:  111111111 

Latches :  1111111111111111110111111111111111110101 
Output :  0 
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4 

Inputs:  111111111 

Latches :  1111111111111111110111111111111111111101 
Output :  0 

5 

Inputs:  111111111 

Latches :  1111111111111111111111111111111111111111 
Output :  1 

Abstract  counterexample: 

000000 

100000 

110001 

111001 

111101 

mill 

Abstract  counterexample  is  invalid  at  step  2 . 

ILP  is: 

94 

Splitting  on  literal  94 

Checking  abstract  model,  hidden  literals  are:  20  22  24  26  28  30  32 
34  36  38  40  42  44  46  48  50  54  56  58  60  62  64  66  68  70  72  74  76  78  80 
84  86  90 

Making  transition  relation 
Done  making  transition  relation 
Preparing  for  model-checking. 

Starting  fixed  point  computation. 

Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Goal  reached 
Goal  Reached  in  5  steps 
Counterexample : 

0 

Inputs:  111111111 

Latches:  0000000000000000000000000000000000000000 
Output :  0 

1 
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Inputs:  111111111 

Latches :  1111111111111111110111111111111011010000 
Output :  0 

2 

Inputs:  111111111 

Latches :  1111111111111111110111111111111111010100 
Output :  0 

3 

Inputs:  111111111 

Latches :  1111111111111111110111111111111111110101 
Output :  0 

4 

Inputs:  111111111 

Latches :  1111111111111111110111111111111111111101 
Output :  0 

5 

Inputs:  111111111 

Latches :  1111111111111111111111111111111111111111 
Output :  1 

Abstract  counterexample: 

0000000 

1000000 

1100100 

1110101 

1111101 

1111111 

Abstract  counterexample  is  invalid  at  step  2 . 

ILP  is: 

90 

Splitting  on  literal  90 

Checking  abstract  model,  hidden  literals  are:  20  22  24  26  28  30  32 
34  36  38  40  42  44  46  48  50  54  56  58  60  62  64  66  68  70  72  74  76  78  80 
84  86 

Making  transition  relation 
Done  making  transition  relation 
Preparing  for  model-checking. 

Starting  fixed  point  computation. 

Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
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Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Goal  reached 
Goal  Reached  in  5  steps 
Counterexample : 

0 

Inputs:  111111111 

Latches:  0000000000000000000000000000000000000000 
Output :  0 

1 

Inputs:  111111111 

Latches :  1111111111111111110111111111111011000000 
Output :  0 

2 

Inputs:  111111111 

Latches :  1111111111111111110111111111111111010000 
Output :  0 

3 

Inputs:  111111111 

Latches :  1111111111111111110111111111111111110100 
Output :  0 

4 

Inputs:  111111111 

Latches :  1111111111111111110111111111111111111101 
Output :  0 

5 

Inputs:  111111111 

Latches :  1111111111111111111111111111111111111111 
Output :  1 

Abstract  counterexample: 

00000000 

10000000 

11010000 

11110100 

11111101 

11111111 

Abstract  counterexample  is  invalid  at  step  3. 

ILP  is: 

84 

Splitting  on  literal  84 

This  refinement  proeess  eontinues,  the  eounterexample  still  being  trivial,  and  the  refinement 
introdueing  more  of  the  latehes  in  the  three-lateh  left  and  right  ehains. 

Quiek  summary:  these  refinements  eliminate  the  “just  turn  the  right  on  and  the  left  on”  abstraet 
eounterexample.  It  takes  multiple  refinement  steps  to  introduee  all  of  the  latehes  in  the  “rightl” 
...  “rights”  and  “leftl”...”left3”  delay  ehains. 
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Goal  Reached  in  5  steps 
Counterexample : 

0 

Inputs:  111111111 

Latches:  0000000000000000000000000000000000000000 
Output :  0 

1 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIOIIIIIIIIIIIIOIIOOOOOO 
Output :  0 

2 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIOIIIIIIIIIIIIIIIOIOOOO 
Output :  0 

3 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIOIIIIIIIIIIIIIIIIIOIOO 
Output :  0 

4 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIOIIIIIIIIIIIIIIIIIIIOI 
Output :  0 

5 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII 
Output :  I 

Abstract  counterexample: 

000000000 

101000000 

111010000 

111110100 

111111101 

111111111 

Abstract  counterexample  is  invalid  at  step  2 . 

ILP  is: 

70  72  74 
70  72  76 
70  74  76 
72  74  76 
70  72  78 
70  74  78 
72  74  78 
70  76  78 
72  76  78 
74  76  78 

No  common  element  in  ILP.  Enter  literals  to  split  on. 
Calling  print_ilp 


Minimized  ILP: 

70  72  74 
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70  72  76 
70  74  76 
72  74  76 
70  72  78 
70  74  78 
72  74  78 
70  76  78 
72  76  78 
74  76  78 


Returned  from  print_ilp 
70  72  74 


At  this  point,  we  have  an  abstract  counterexample  that  effectively  says  “turn  the  right  controller 

o 

on,  turn  the  left  controller  on,  wait  for  three  cycles,  then  turn  the  valid  latch  on.”  Now  we  get 
our  first  interesting  ILP,  as  the  solver  generates  a  spurious  counterexample,  the  sensing  logic 
enabling  the  left  controller  to  turn  on  is  still  abstract.  In  fact,  the  controller  can  only  be  turned  on 
if  the  sensors  detect  the  right  condition  and  enough  of  the  redundant  sensors  vote  to  enable  it. 
There  are  a  large  number  of  combinations  of  possible  sensor  result  latches  that  could  lead  to  this 
voting  outcome,  and  the  D-Wave  ILP  solver  chooses  a  minimal  covering  set  (literals  70,  72,  74). 


Goal  reached 

Goal  Reached  in  6  steps 

Counterexample : 

0 

Inputs:  111111111 

Latches:  0000000000000000000000000000000000000000 
Output :  0 

1 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIOIIIIIIOOOIIIOIIOOOOOO 
Output :  0 

2 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIOIIIIIIIIIIIIIOIOIOOOO 
Output :  0 

3 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIOIIIIIIIIIIIIIIIIOOIOO 
Output :  0 

4 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIOIIIIIIIIIIIIIIIIIIOOI 


g 

Since  the  logic  behind  computing  validity  is  still  abstract,  the  system  considers  that  it  could  just 
turn  the  valid  latch  on. 

APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 

46 


5 


Output :  0 


Inputs:  111111111 

Latches :  1111111111111111110111111111111111111110 
Output :  0 

6 

Inputs:  111111111 

Latches :  1111111111111111111111111111111111111111 
Output :  1 

Abstract  counterexample: 

000000000000 

100001000000 

111110010000 

111111100100 

111111111001 

111111111110 

111111111111 

Abstract  counterexample  is  invalid  at  step  4 . 

ILP  is: 

36  38  40  56 
40  56 
36  38  56 
56 

Splitting  on  literal  56 

Checking  abstract  model,  hidden  literals  are:  20  22  24  26  28  30  32 
34  36  38  40  42  44  46  48  50  54  58  60  62  64  66  68  76  78  80  86 
Making  transition  relation 
Done  making  transition  relation 
Preparing  for  model-checking. 

Starting  fixed  point  computation. 

Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Goal  reached 
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Goal  Reached  in  6  steps 
Counterexample : 

0 

Inputs:  111111111 

Latches:  0000000000000000000000000000000000000000 
Output :  0 

1 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIOOIIIIIIOOOIIIOIIOOOOOO 
Output :  0 

2 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIOOIIIIIIIIIIIIIOIOIOOOO 
Output :  0 

3 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIOOIIIIIIIIIIIIIIIIOOIOO 
Output :  0 

4 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIOOIIIIIIIIIIIIIIIIIIOOI 
Output :  0 

5 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIOIIIIIIIIIIIIIIIIIIIIO 
Output :  0 

6 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII 
Output :  I 

Abstract  counterexample: 

0000000000000 

1000001000000 

1011110010000 

1011111100100 

1011111111001 

1011111111110 

1111111111111 

Abstract  counterexample  is  invalid  at  step  3. 

ILP  is: 

36  38  40  54 
40  54 
36  38  54 
54 

Splitting  on  literal  54 

The  next  eounterexample  has  the  left  sensing  on,  and  a  eoneurrent  enabling  of  the  right  channel. 
The  model-checker  detects  that  this  abstract  counterexample  is  incorrect,  because  of  the  trace 
validity  constraints  that  capture  the  mutual  exclusion  between  the  right  and  left  controller 
APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 

48 


enabling  conditions.  In  order  to  check  this  correctly,  the  model  must  be  refined  to  take  into 
account  a  combination  of  different  bits  of  validity  logic,  including  the  one-shot  logic  that 
maintains  the  validity  of  the  trace,^  and  the  left  integrity  constraints.  The  next  two  steps  in  this 
refinement  are  trivial:  we  add  the  one-shot  validity  logic  to  the  model. 

This  step  adds  the  second  latch  required  for  the  one-shot  validity  logic. 

Goal  Reached  in  6  steps 
Counterexample : 

0 

Inputs:  111111111 

Latches:  0000000000000000000000000000000000000000 
Output :  0 

1 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIOOIOOIIIIIIOOOIIIOIIOOOOOO 
Output :  0 

2 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIOOIOOIIIIIIIIIIIIIOIOIOOOO 
Output :  0 

3 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIOOIOOIIIIIIIIIIIIIIIIOOIOO 
Output :  0 

4 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIOOIIIIIIIIIIIIIIIIIIOOI 
Output :  0 

5 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIOIIIIIIIIIIIIIIIIIIIIO 
Output :  0 

6 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII 
Output :  I 

Abstract  counterexample: 

00000000000000 

10000001000000 

10011110010000 

10011111100100 

10011111111001 

11011111111110 

11111111111111 

Abstract  counterexample  is  invalid  at  step  2 . 

ILP  is: 


^  If  an  invalid  state  is  detected  in  the  trace,  then  it  is  always  invalid. 
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36  38  40  48 
36  38  42  48 
36  38  48 
40  44  48 
40  48 
42  44  48 
44  48 
42  48 
48 

36  38  40  50 
36  38  40 
36  38  42  50 
36  38  50 
36  38  42 
40  44  50 
40  50 
40  44 
42  44  50 
44  50 
42  50 
50 

42  44 

No  common  element  in  ILP.  Enter  literals  to  split  on. 
Calling  print_ilp 


Minimized  ILP: 
48 

36  38  40 
36  38  42 
40  44 
50 

42  44 


Returned  from  print_ilp 
36  38  44  48  50 

At  this  point,  the  model  has  been  refined  to  include  enough  variables  to  recognize  that  the 
validity  flag  starts  on,  and  if  turned  off  cannot  be  restored.  This  eliminates  spurious 
counterexamples  where  the  adversary  takes  the  system  through  invalid  transitions  to  turn  on  the 
right  and  left  controllers  simultaneously,  and  then  flips  the  valid  flag  on  from  off  The  new 
abstract  counterexamples  involve  assignments  to  the  sensor  inputs  that  would  allow  the  aircraft 
to  enable  both  right  and  left  controllers  simultaneously.  There  are  no  valid  assignments  that  do 
so,  and  the  model  checker  detects  that  the  counterexamples  are  spurious.  The  final  ILP  requires 
the  solver  to  identify  a  covering  set  of  variables  that  is  sufficient  to  capture  enough  of  the 
validity-checking  logic  to  eliminate  all  of  these  abstract  counterexamples.  It  chooses  enough  of 
the  variables  to  enforce  a  valid  assignment  of  sensor  variables  to  the  left  side  logic,  and  enough 
to  enforce  the  consistency  relations  between  the  left  and  right  sides. 
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Goal  reached 

Goal  Reached  in  6  steps 

Counterexample : 

0 

Inputs:  110111111 

Latches:  0000000000000000000000000000000000000000 
Output :  0 

1 

Inputs:  IIOIIIIII 

Latches :  IIIIIIIIOOIOOOOIIOOIIIIIIOOOIIIOIIOOOOOO 
Output :  0 

2 

Inputs:  IIOIIIIII 

Latches :  IIIIIIIIOOIOOOOIIOOIIIIIIOOIIIIIOIOIOOOO 
Output :  0 

3 

Inputs:  IIIIIIIOI 

Latches :  IIIIIIIIOOIOOOOIIOOIIIIIIOOIIIIIIIIOOIOO 
Output :  0 

4 

Inputs:  IIIIIIIOI 

Latches :  IIIIIIIIIIIIIIIOIOOIIIIIIOOIIIIIIIIIIOOI 
Output :  0 

5 

Inputs:  IIIIIIIOI 

Latches :  IIIIIIIIIIIIIIIOIIOIIIIIIIIIIIIIIIIIIIIO 
Output :  0 

6 

Inputs:  IIIIIIIII 

Latches :  IIIIIIIIIIIIIIIOIIIIIIIIIIIIIIIIIIIIIIII 
Output :  I 

Abstract  counterexample: 

0000000000000000000 

0000110000001000000 

0000110000110010000 

0000110000111100100 

1111010000111111001 

1111011011111111110 

1111011111111111111 

Abstract  counterexample  is  invalid  at  step  2 . 

ILP  is: 

40 

42 

No  common  element  in  ILP.  Enter  literals  to  split  on. 
Calling  print_ilp 


Minimized  ILP: 

40 

42 
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Returned  from  print_ilp 
40  42 

Will  split  on: 

40  42 
OK?  y 

Splitting  on  40  42 

Checking  abstract  model,  hidden  literals  are:  20  22  24  26  28  30  32 

34  46  58  60  62  64  66  68  76  78  80  86 

Making  transition  relation 

Done  making  transition  relation 

Preparing  for  model-checking. 

Starting  fixed  point  computation. 

Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Checking  for  goal 
Computing  successors 
Checking  for  new  state 
Reached  fixpoint  after  7  steps 
Goal  not  reached 

No  abstract  path  to  goal  state:  circuit  is  safe. 


This  is  almost  enough  to  eliminate  all  of  the  spurious  counterexamples.  The  final  refinement 
requires  the  addition  of  two  more  variables  (40  and  42),  which  round  out  some  of  the  input 
consistency  constraints,  and  then  the  solver  establishes  that  the  failure  state  is  unreachable,  and 
the  process  concludes. 


4,6  Evidence  for  quantum  behavior  in  the  DW2  processor 

Even  though  it  was  not  an  explicit  part  of  the  QCHECK  proposed  research,  we  would  like  to 
present  very  important  results  regarding  the  quantum  nature  of  the  D-Wave  processor  that  were 
obtained  through  research  conducted  in  parallel  to  the  main  QCHECK  work. 
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The  D-Wave  processor  was  designed  to  exploit  quantum  mechanical  effects  in  solving 
combinatorial  optimization  problems.  However,  in  order  to  make  the  design  scalable,  certain 
compromises  in  the  hardware  had  to  be  made.  The  main  result  of  these  was  that  only  a  restricted 
set  of  measurements  were  allowed  to  investigate  the  state  of  the  processor,  and  hence  we  lack  an 
important  set  of  tools  to  determine  if  the  state  of  the  system  indeed  behaves  quantum 
mechanically.  We  worked  on  developing  two  approaches  to  address  this  issue:  in  one,  we 
designed  a  very  simple  problem  whose  output  statistics  were  different  for  a  simple  classical 
model  and  a  quantum  model;  in  the  other,  using  partial  state  information  about  the  state  of  the 
processor  during  the  annealing  we  proved,  using  a  novel  analytical  technique  developed  for  this 
purpose,  that  the  state  of  the  system  has  entanglement  at  some  points  during  the  annealing.  We 
give  a  brief  description  of  both  approaches. 

4.6.1  Quantum  signature 

We  constructed  a  simple  toy  problem  involving  8  qubits  that  had  the  following  characteristics:  (i) 
the  ground  state  was  IV^’^-fold  degenerate;  (ii)  the  ground  states  were  divided  in  two  sets,  a 
cluster  whose  16  states  could  be  transformed  to  another  ground  state  by  flipping  a  single  qubit, 
and  another  single  state  that  required  flipping  4  qubits  in  order  to  obtain  another  ground  state; 
(iii)  the  energy  landscape  of  the  problem  had  no  local  minima.  We  studied  the  evolution  of  the 
populations  of  the  17  ground  states  according  to  a  simulated  annealing  evolution  (which  is  a 
simple  classical  model  of  a  thermalizing  system)  and  found  that  the  isolated  state  initially  had  a 
higher  population  than  the  states  in  the  cluster.  On  the  other  hand,  a  quantum  mechanical  model 
of  the  system  predicted  that  the  isolated  state  would  have  a  suppressed  population  with  respect  to 
the  cluster  states.  We  then  performed  the  experiment  on  the  Rainier  processor  (a  128-qubit 
processor,  that  preceded  DW2),  and  found  that  the  experimental  statistics  agreed  with  the 
quantum  model  and  disagreed  with  the  classical  model.  These  results  were  published  in  Nature 
Communications  [18]  (see  appendix  A). 

4.6.2  Evidence  of  entanglement 

The  second  approach  was  a  collaboration  with  D-Wave  to  analyze  the  entanglement  of  a  set  of  8 
qubits  in  the  DW2  processor  (Vesuvius  chip).  DW2  has  extra  controls  that  allow  for  certain 
measurements  that  were  not  available  for  the  Rainier  processor  (DWl).  Using  these  extra 
controls,  an  experiment  was  designed  in  which  one  qubit  was  used  as  a  probe  to  measure  some 
elements  of  the  state  of  a  set  of  eight  interacting  qubits.  The  extra  controls  allowed  for 
independent  annealing  of  the  probe  and  the  system  under  study.  With  this  setup,  researchers  from 
D-Wave  measured  the  populations  of  the  ground  state  and  the  first  excited  state  of  the  8-qubit 
system  at  different  points  during  the  annealing.  These  two  quantities  are  just  2  of  the  65535  real 
parameters  needed  to  completely  determine  the  quantum  state  of  the  system.  We  at  Information 
Sciences  Institute  (ISI)  developed  a  new  technique  to  analyze  this  information  that  allowed  us  to 
conclude  that  no  matter  what  the  values  of  the  unknown  parameters  were,  any  state  that  was 
compatible  with  the  measured  values  for  the  two  populations  had  to  be  entangled.  This  was  the 
first  experimental  proof  of  the  presence  of  entanglement  in  the  D-Wave  processor,  and  it 
provides  definite  evidence  of  the  quantum  mechanical  nature  of  the  device.  This  work  was 
accepted  for  publication  in  Physical  Review  X,  but  has  not  been  published  at  the  time  this  report 
is  being  written.  A  preprint  [19]  is  attached  in  appendix  A. 
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5  CONCLUSIONS 

The  goal  of  the  QCHECK  project  was  to  study  the  feasibility  of  integrating  the  capabilities  of  the 
D-Wave  adiabatic  quantum  annealer  into  a  model  checking  framework  for  the  certification  of 
complex  systems.  The  results  of  this  project  show  that  this  is  indeed  possible.  By  identifying  a 
particular  approach  to  model  checking  that  required  solving  combinatorial  optimization 
problems,  we  were  able  to  off  load  these  hard  computational  tasks  to  the  DW2  processor, 
construct  a  hybrid  algorithm  that  applied  classical  model  checking  techniques,  and  call  upon  the 
DW2  processor  whenever  needed.  The  CEGAR  approach  to  model  checking  was  a  very  good  fit 
to  exploit  the  strengths  of  the  DW2  device.  In  fact,  the  main  difficulty  encountered  in  this  part  of 
the  project,  was  the  lack  of  a  model  checking  package  that  was  open  enough  to  allow  us  to 
extract  the  required  combinatorial  optimization  problems  that  needed  to  be  solved  by  DW2. 
These  packages  typically  use  their  own  solvers,  and  do  not  provide  the  user  with  the  low  level 
access  to  extract  them.  But  looking  into  the  future,  this  project  shows  that  the  DW2  processor 
can  be  trivially  integrated  into  the  CEGAR  approach,  if  the  necessary  access  is  provided. 
Eurthermore,  any  advances  in  the  model  checking  side  and  DW2  side  will  not  alter  the  basics  of 
this  integration,  so  improvements  and  new  developments  in  both  components  will  only  improve 
the  integrated  approach. 

One  of  the  main  reasons  for  our  interest  in  combining  the  DW2  device  with  a  model  checking 
approach  is  the  possibility  that  computational  speedups  may  be  provided  by  the  adiabatic 
quantum  optimization  approach.  The  question  of  quantum  speedups  is  not  easy  to  answer  at  this 
point  in  time.  In  particular,  the  size  of  the  devices  currently  available  may  not  be  large  enough  to 
fully  break  beyond  the  capabilities  of  classical  computers.  Eurthermore,  it  is  very  difficult  to 
prove  that  a  quantum  device  is  better  than  all  possible  classical  devices  at  solving  a  particular 
class  of  problems.  We  are  then  left  with  comparing  the  quantum  device  against  a  set  of  classical 
algorithms  on  a  particular  ensemble  of  problems,  and  it  is  not  clear  if  this  performance  will  carry 
to  more  general  problems,  or  if  another  more  efficient  classical  algorithm  will  be  developed.  In 
this  project  we  compared  the  performance  of  DW2  against  a  heuristic  classical  solver  that  has 
performed  well  in  several  algorithmic  competitions.  Even  though  DW2  performed  well  against 
it,  it  is  not  clear  what  conclusions  we  can  draw  at  this  point,  since  the  classical  solver  was 
designed  for  general  problems  and  not  optimized  to  the  architecture  of  DW2.  The  question  of 
speedup  has  been  (and  still  is)  extensively  studied  by  many  groups  around  the  world,  and  the 
results  are  still  inconclusive  [20].  On  the  other  hand,  the  results  are  also  promising  in  the  sense 
that  the  performance  of  DW2  is  on  par  with  very  optimized  algorithms  running  on  very  fast 
hardware  [20]. 

In  this  project  we  developed  some  of  the  tools  that  would  be  necessary  for  any  practical 
application  involving  DW2.  One  of  the  main  obstacles  in  using  the  processor  is  the  limited 
connectivity  of  the  underlying  architecture  that  severely  limits  the  type  of  problems  that  can  be 
embedded  exactly.  It  is  thus  unavoidable  that  users  will  have  to  rely  on  approximate  or  heuristic 
embedding  techniques  to  solve  a  broader  set  of  problems.  This  has  its  own  set  of  drawbacks  as 
the  solutions  we  found  are  only  approximate,  and  it  is  not  clear  a  priori  how  good  those 
approximations  are.  The  sequential  embedding  technique  developed  in  this  project  shows  one 
possible  way  forward  and  the  results  are  encouraging  for  optimization  problems,  but  not  so  much 
for  decision  problems.  The  issue  is  that  decision  problems  are  implemented  in  DW2  by  first 
encoding  them  as  optimization  problems:  the  answer  to  the  decision  problem  is  YES  if  the 
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minimum  of  a  certain  discrete  function  is  equal  to  a  certain  value  (usually  zero),  and  NO 
otherwise.  The  problem  with  this  type  of  formulation  is  that  the  DW2  processor  does  not  provide 
a  guarantee  that  the  ground  state  has  been  found.  So,  if  the  best  configuration  provided  by  the 
quantum  processor  gives  us  the  known  minimum  of  the  function,  we  know  the  answer  to  the 
decision  problems  is  YES,  but  if  it  is  not  equal  to  that  minimum,  we  cannot  say  that  the  answer  is 
NO,  because  there  is  no  guarantee  that  the  answer  found  is  the  best  possible  one.  We  faced  this 
issue  when  attempting  to  use  DW2  to  check  the  validity  of  abstract  counterexamples  (i.e., 
checking  is  there  was  a  corresponding  counterexample  in  the  original  system).  This  required 
solving  a  SAT  decision  problem,  and  the  fact  that  we  could  not  guarantee  that  the  best  solution 
was  found  would  lead  us  to  a  wrong  assertion  about  our  problem  (i.e.,  that  a  real  counter 
example  existed  when  that  was  not  true). 

Another  important  feature  that  was  shown  in  parallel  with  the  work  performed  for  QCHECK  was 
the  quantum  mechanical  behavior  of  the  D-Wave  processor.  We  developed  two  approaches,  one 
that  gave  evidence  of  a  quantum  signature  (but  fell  short  of  being  a  proof),  and  another  one  that 
provided  experimental  evidence  of  the  existence  of  entangled  state  during  the  annealing  process. 
This  last  one  constitutes  a  definite  proof  of  the  quantum  mechanical  nature  of  the  device  and  it  is 
a  very  important  step  forward  in  the  field  of  practical  implementations  of  programmable 
quantum  devices.  There  is,  however,  a  lot  of  work  to  be  done  to  determine  if  these  quantum 
effects  are  being  exploited  in  a  way  that  provides  a  computational  speedup. 

In  summary,  integrating  adiabatic  quantum  optimization  techniques  with  a  model  checking 
approach  is  feasible  and  relatively  straightforward.  The  framework  developed  in  this  project 
could  be  easily  followed  in  the  future  provided  the  model  checking  packages  used  provide  the 
user  with  the  required  information.  Developing  newer  techniques  to  go  beyond  the  connectivity 
limitations  of  the  DW2  device  will  certainly  help  improve  the  performance  of  the  approach.  And 
hopefully,  newer  generations  of  the  device  will  provide  the  desired  computational  speedup  and 
transform  this  proof  of  concept  into  a  practical  tool. 
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7  APPENDIX  A  -  Publications  and  Presentations 


The  following  two  publications,  regarding  the  quantum  nature  of  the  D-Wave  device,  were 
written  in  parallel  with  the  work  performed  for  the  QCHECK  project.  Although  not  a  part  of  the 
proposed  research,  these  results  are  included  to  complement  and  give  better  context  to  the  results 
presented  in  this  report. 


•  Quantum  Signature 

S.  Boixo,  T.  Albash,  F.M.  Spedalieri,  N.  Chancellor,  D.  A.  Lidar;  “Experimental  signature  of 
programmable  quantum  annealing”;  Nature  Communications  4,  Article  number:  2067, 
published  online  28  June  2013. 


•  Quantum  entanglement 

T.  Canting,  A.J.  Przybysz,  A.  Yu.  Smirnov,  F.M.  Spedalieri,  M.H.  Amin,  A.J.  Berkley,  R. 
Harris,  F.  Altomare,  S.  Boixo,  P.  Bunyk,  N.  Dickson,  C.  Enderud,  J.P.  Hilton,  E.  Hoskinson, 
M.W.  Johnson,  E.  Eadizinsky,  N.  Eadizinsky,  R.  Neufeld,  T.  Oh,  I.  Perminov,  C.  Rich,  M.C. 
Thom,  E.  Tolkacheva,  S.  Uchaikin,  A.B.  Wilson,  G.  Rose;  “Entanglement  in  a  quantum 
annealing  processor”;  arXiv:  1401. 3500,  accepted  for  publication  in  Physical  Review  X, 
2014. 
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8  APPENDIX  B  -  Description  of  CEGAR-DW2  integration  code 


In  this  appendix  we  provide  a  brief  deseription  of  the  main  eode  eomponents  needed  to  integrate 
the  CEGAR  approaeh  with  the  DW2  eapabilities.  As  diseussed  in  Seetion  4.4,  we  were  foreed  to 
eonstruet  a  simple  model  eheeking  program  (based  on  freely  available  BDD  libraries)  beeause 
publiely  available  model  eheeking  paekages  did  not  provide  the  user  with  low  level  information 
needed  to  eonstruet  the  ILPs  required  during  the  CEGAR  proeess.  This  eode  (and  supporting 
fdes)  ean  be  aeeessed  at  http://www.isi.edU/people//fspedali/QCHECK_eode 

To  run  the  integrated  CEGAR  proeess,  we  eall  a  Perl  seript  named  cegar-solver- 
loop.pl.  This  seript  takes  as  an  argument  the  name  of  the  file  where  the  system  is  deseribed  in 
AIG  format.  The  CEGAR  proeess  is  then  applied:  an  abstraetion  is  generated,  the  model  eheeker 
is  ealled  to  eheek  the  required  property,  and  if  a  eounterexample  is  found,  it  is  eheeked  using  a 
regular  SAT  solver.  When  the  eounter  example  is  found  to  be  spurious,  an  ILP  is  generated  to 
solve  the  Minimal  Separating  Set  problem  that  would  allow  us  to  refine  the  abstraetion.  Onee  the 
lEP  is  generated,  it  is  written  to  a  temporary  file.  Then  the  DW2  solver  is  ealled  (see  the  line 

my  $answer  ='matlab  -nodisplay  -nosplash  -r 
"SolverILP ( ' $tmpf lie ' ) "  |  tail  -n  +11'; 

in  cegar-solver-loop.pl).  This  eommand  ealls  Matlab  and  runs  the  funetion 
SolverILP,  that  takes  the  name  of  the  file  where  the  ILP  is  stored  as  an  argument,  and  then 
ealls  other  funetions  that  write  the  ILP  as  a  QUBO  problem  and  solves  it  using  the  heuristie 
embedding  eoded  in  the  funetion  SEBREMforQUBO.  The  answer  is  sent  to  the  standard  output, 
where  the  Perl  seript  grabs  it  and  eontinues  with  the  CEGAR  proeess,  refining  the  abstraetion 
aeeording  with  the  solution  of  the  ILP.  The  eode  ean  be  run  in  verbose  mode,  where  the 
sueeessive  ILPs  solved  are  shown.  A  detailed  transeript  of  one  sueh  run  was  presented  in  Seetion 
4.5. 
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9  LIST  OF  SYMBOLS,  ABBREVIATIONS,  AND  ACRONYMS 


G 

Pauli  operator 

2-SAT 

2-Satisfiability  problem 

AIG 

And  Inverter  Graph 

AIGER 

Eormat  for  And-Inverter  Graphs 

API 

Applieation  Programming  Interfaee 

AX-MAXSAT 

Exaet  MAX-SAT  solver 

BDD 

Binary  Deeision  Diagram 

BMC 

Bounded  Model  Cheeking 

CE 

Counter  Example 

CEGAR 

Counter-example  Guided  Abstraetion  Refinement 

CMU 

Carnegie  Mellon  University 

CNE 

Conjunetive  Normal  Eorm 

CTE 

Computational  Tree  Eogie 

CUDD 

CU  Deeision  Diagram  paekage 

DIMACS 

Computer-readable  format  for  Boolean  satisfiability  problems 

DWl 

D-Wave  One 

DW2 

D-Wave  Two 

EBK 

Eondazione  Bruno  Kessler 

ECS 

Elight  Control  System 

h 

Veetor  of  loeal  magnetie  fields 

lEP 

Integer  Einear  Program 

ISI 

Information  Scienees  Institute 

J 

Matrix  of  inter  qubit  eouplings 

laig 

Eisp-fiavored  AIG 

ETE 

Einear  Temporal  Eogie 

MAX-2-SAT 

Maximum  2-Satisfiability  problem 

MaxWalkSAT 

Walks  AT  variant  for  weighted  SAT  problem  solver 

MC 

Model  Checking 
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NP 

Non-deterministic  Polynomial  class  of  decision  problems 

NP-hard 

Computational  complexity  class  that  contains  the  elass  NP 

NuSMV 

New  Symbolie  Model  Checker  for  System  Verification 

NuXMV 

New  Symbolic  Model  Cheeker  for  the  Analysis  of  Synehronous 
finite-state  and  infinite-state  Systems 

OBDD 

Ordered  Reduced  Binary  Deeision  Diagram  paekage 

QCHECK 

A  Quantum  Computing  Approach  to  Model  Checking  for 
Advanced  Manufaeturing  Problems 

QUBO 

Quadratic  Unconstrained  Binary  Optimization 

rf-SQUID 

Radio  Erequency  Superconducting  Quantum  Interference  Device 

SAT 

Boolean  satisfiability 

SEBREM 

Sequential  Embedding  By  Relative  Entropy  Minimization 

SMV 

Symbolic  Model  Cheeker  for  System  Verification 

SQUID 

Superconducting  Quantum  Interference  Deviee 

WalkSAT 

Boolean  satisfiability  problem  solver 

XME 

Extended  Markup  Eanguage 
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